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Summary 

A  fuel  cell  is  an  electrochemical  energy  conversion  device  that  combines  hydrogen  and  oxygen  (from  the 
air)  to  produce  electricity.  Currently,  the  Proton  Exchange  Membrane  technology  (PEM)  is  the  most 
promising  and  is  being  developed  by  United  Technologies  Fuel  Cells.  One  of  the  major  obstacles  in  the 
commercialization  of  PEMFC  technology  is  the  performance  loss  at  high  current  density  due  to  liquid 
water  accumulation  in  the  porous  components  of  the  cathode  compartment.  This  accumulation  limits 
oxygen  transport  to  the  cathode  and  reduces  PEMFC  efficiency  and  performance  stability. 

In  this  work  we  aim  developing  a  novel  theoretical  framework  to  evaluate  the  feasibility  of  attaining 
significant  improvement  of  fuel  cells  performance  and  stability  by  increasing  the  transport  processes  in 
porous  partially  fluid  filled  cathode  compartment  through  the  application  of  acoustic  and  structural 
waves.  The  specific  goal  of  the  project  is  to  develop  a  generic  model  of  partially  saturated  multiphase 
flow  in  highly  porous  media  (Fuel  Cell  cathode  compartment)  under  structural  waves. 

We  have  developed  a  unified  model  of  structural/acoustic  wave  propagation  in  the  PEM  cathode 
compartment  and  coupled  it  with  mass  transfer  in  the  porous  media.  This  is  the  first  known  to  authors 
model  of  such  kind  fully  coupling  acoustics  and  filtration  in  porous  medium.  A  new  perturbation 
technique  allows  deriving  the  governing  equations  for  multiphase  flow  under  acoustic  streaming.  It  has 
been  demonstrated  that  the  phase  saturations  have  strong  impact  on  wave  dynamics  in  porous  continuum. 
Explicit  expressions  for  generalized  multiphase  Biot-type  coefficients  that  explicitly  depend  on  the  phase 
saturations  have  been  obtained.  A  novel  generalized  filtration  law  that  accounts  for  dynamic  loadings, 
varying  saturation,  and  solid  structure  distortion  describing  mass  transfer  in  this  complex  but  generic 
system  has  been  found.  We  also  have  been  able  to  connect  hysteretic  phenomena  to  relaxation  process 
due  to  filtration  velocity  gradient.  The  Lagrangian  approach  has  been  employed  in  order  to  develop 
consistent  government  equations.  We  have  also  conducted  highly  precise  tests  on  flow  measurements 
throughout  porous  media  with  and  without  applied  acoustic  signals. 

The  developed  homogenization  technique  allows  deducing  effective  body  forces  in  mass  transfer 
equations  from  acoustic  terms  and  simple  analytical  evaluation  of  frequency  and  wave  number 
dependencies  of  the  effective  mass  driving  body  forces.  It  has  been  demonstrated  using  averaging 
principle  that  vibration  gives  rise  to  net  change  of  saturation  inside  porous  medium.  This  effect  is 
analogous  to  the  known  phenomenon  of  acoustic  streaming  but  has  not  been  discovered  before. 

Numerical  results  demonstrate  that  applied  structural  vibration/acoustic  loading  (i)  intensifies  the  process 
dynamics  or  in  other  words  decreases  time  of  transient  process;  (n)  might  drain  out  up  to  5-10%  of 
distributed  in  porous  media  water  if  the  mechanism  preventing  reverse  acoustic-driven  imbibitions  is  in 
place;  {Hi)  The  fast  water  removal  from  the  surface  is  crucial  for  acoustic-driven  drying.  Based  on  the 
numerical  and  experimental  results  number  of  practical  recommendations  optimizing  material  selections 
and  performance  regime  has  been  made.  Developed  methodology  is  useful  for  wide  range  of  Fuel  Cell 
problems  as  well  as  for  wide  range  of  other  porous  structures  and  could  serve  as  an  important  design  tool. 
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Program  Overview 


Technical  Objectives 

In  this  work  our  target  is  to  develop  a  theoretical  framework  to  evaluate  the  feasibility  of 
attaining  significant  improvement  of  fuel  cells  performance  and  stability  by  increasing  the 
transport  processes  in  porous  fluid  filled  cathode  compartment  through  the  application  of 
acoustic  and  structural  waves.  The  specific  goal  of  the  project  is  to  develop  a  generic  model  of 
unsaturated  multiphase  flow  in  highly  porous  media  (Fuel  Cell  cathode  compartment)  under 
structural  waves. 

In  order  to  attain  this  goal  the  novel  model  predicting  effects  of  structural  acoustics  on  partially 
saturated  mass  transfer  as  well  as  effect  of  dynamic  saturation  on  wave  characteristics  should  be 
developed. 

To  the  best  of  our  knowledge,  structural  or  acoustic  excitation  has  never  been  applied  to  the  fuel 
cells  in  particular  and  to  the  porous  membranes  in  general.  This  project  combines  computational 
and  experimental  efforts  including: 

1)  Developing  of  fully  coupled  physics-based  model  of  structural  and  acoustic  waves 
propagation  in  partially  water-  filled  porous  media  and  induced  multiphase  (gas/water) 
transport  through  the  partially  water-  filled  vibrating  porous  media; 

2)  Experimental  calibration  and  verification  of  the  models  using  company-developed  rigs 
and  fuel  cells; 

3)  Numerical  demonstration  of  transport  process  and  fuel  cell  performance  enhancement. 

Technical  Approach 

We  aim  to  improve  the  efficiency  of  a  PEM  fuel  cell  by  intensifying  the  gas/water  transport 
processes  using  either  structural  or  acoustic  waves.  Structural  excitation  can  be  used  for  inducing 
waves  that  would  propagate  within  the  different  layers  of  the  fuel  cell  and  perturb  the  trapped 
liquid  water  within  pores,  thereby  enhaneing  the  transport  process.  Introducing  pulsations  in  the 
airflow  past  the  cathode  can  also  intensify  the  diffusion  process;  this  is  expected  to  set  up 
pressure  gradients  in  the  fuel  cell  that  would  pull  the  trapped  water  cavities  into  the  air  stream. 

The  first  step  in  this  effort  involves  choosing  of  the  modeling  strategy.  All  possible  approaches 
can  be  separated  into  two  categories:  1.  Continuum  mechanics  analysis  including  wave 
propagation  and  mass  transfer  analysis;  and  2)  Micro-structural  /  micro-mechanical 
considerations  accounted  for  the  local  interactions.  The  second  approach  is  good  for  the  analysis 
of  the  role  of  material  morphology  but  not  suitable  for  waves  and  flux  analysis.  We  chose  the 
macroscopic  (continuum  mechanics)  modeling  approach,  which  is  the  best  suitable  for  the 
analysis  of  wave  propagations,  multiphase  fluxes  (not  one  or  even  hundreds  droplets)  in  real 
geometry  and  for  coupling  acoustic  and  mass  transfer  properties  and  time-scales  The 
macroscopic  modeling  framework  naturally  (i)  provides  a  unified  system  of  governing  equations 
describing  the  whole  physical  system;  (ii)  provides  the  solution  of  problems  with  real  geometry. 


7 


layers,  configuration,  regimes,  and  boundary  conditions;  (iii)  investigates  effect  combining  both 
“elastic-wave”  part  and  “viscous  -  mass  transfer”  part  coupling  transport  and  acoustic 
properties;  (iv)  the  dynamic  model  provides  a  comprehensive  description  of  multiphase  flux, 
amount  of  trapped  water  and  system  behavior  during  drainage  and  imbibitions  -  hysteresis. 

The  second  step,  which  will  be  reported,  is  the  development  of  the  novel  governing  equations 
accounted  for  multiphase  unsaturated  flow  in  vibrating  porous  media.  The  system  should  be 
developed  in  a  consistent  way  to  accommodate  simultaneously  all  physical  phenomena.  Using 
the  developed  governing  system  (i)  We  study  wave  propagation  and  attenuation  in  partially 
saturated  porous  media,  (ii)  We  calculate  different  waves  in  the  structure  and  couple  them  with 
pressures  and  water/gas  fluxes,  (iii)  We  calculate  inhomogeneous  water  saturation  level  with 
space  and  time  -  S(x,t)  and  relate  it  with  oxygen  transport,  (iv)  This  approach  needs  to  be 
matched  with  micromechanics  in  order  to  complete  the  governing  system  with  constitutive  laws. 

This  is  followed  by  validation  of  the  models  in  developed  controlled  laboratory  experiments. 
We  conduct  highly  precise  tests  on  flow  measurements  throughout  porous  media  with  and 
without  applied  acoustic  signals.  Also,  we  conduct  the  full  scale  test  on  company  developed  Fuel 
Cell  in  order  to  compare  power  output  with  and  without  acoustic  impact  and  under  different 
acoustic  regimes.  These  tests  serve  as  model  calibration  and  validation. 

The  validated  models  will  then  be  used  for  predicting  the  level  of  structural/acoustic  response 
needed  for  achieving  the  targeted  improvements  in  fuel  cell  performance.  The  high-level 
actuation  methodology  will  be  chosen  based  on  energy  and  efficiency  considerations.  The  choice 
of  actuation  scheme  will  then  be  used  to  identify  promising  actuator  design. 

The  major  accomplishments  are  as  follows: 

(i)  We  have  developed  a  unified  model  of  structural/acoustic  wave  propagation  in  the  PEM 
cathode  compartment  and  coupled  it  with  multiphase  mass  transfer  in  the  partially  saturated 
porous  media.  This  is  the  first  successful  attempt  to  model  coupling  between  acoustics  and 
filtration  in  porous  massive. 

(ii)  It  has  been  demonstrated  that  the  phase  saturations  have  strong  impact  on  wave  dynamics  in 
porous  continuum.  Explicit  expressions  for  generalized  multiphase  Biot-type  coefficients  that 
explicitly  depend  on  the  phase  saturations  have  been  obtained.  We  also  have  been  able  to 
connect  hysteretic  phenomena  to  relaxation  process  due  to  filtration  velocity  gradient. 

(iii)  A  novel  generalized  filtration  law  (generalized  D’Arcy  law)  that  accounts  for  dynamic 
loadings,  varying  saturation,  and  solid  structure  distortion  describing  mass  transfer  in  this 
complex  but  generic  system  has  been  developed. 

(iv)  The  new  developed  homogenization  technique  is  crucial  for  numerical  analysis  of  high 
frequency  impacts.  It  allows  circumventing  the  restriction  on  maximum  allowed  number  of  finite 
elements  to  resolve  spatial  frequency  of  the  acoustic  wave.  Using  asymptotic  analysis  we  obtain 
simple  analytical  evaluations  of  frequency  and  wave  number  dependencies  of  the  effective  mass 
driving  forces 

(v) .  We  conduct  the  extremely  precise  experimental  procedure  to  measure  bi-layer/substrate 
permeability  under  vibration.  Test  results  for  different  substrates  demonstrate  (i)  existence  of  the 
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effect,  (ii)  instability  in  highly  hydrophobic  structures  under  vibration  and  the  need  in  wave 
generation. 

(vi)  Extensive  numerical  analysis  of  the  vibration-induced  mass  transfer  has  been  conducted  and 
on  the  base  of  obtained  numerical  results  (as  well  as  experimental  data)  we  have  made  the 
practical  conclusions  focused  on  (i)  industrial  (fuel  cells)  implementation  and  (ii)  optimal 
material  design  rwequirements  and  operational  conditions  (iii)  fundamental  understanding  of 
mass  transfer  in  unsaturated  porous  media. 

Developed  methodology  is  useful  for  wide  range  of  Fuel  Cell  problems  as  well  as  for  wide  range 
of  other  porous  structures  and  could  serve  as  an  important  design  tool. 


Technical  Accomplishments 

1 .  Technical  Background  and  Problem  Formulation 


Figure  1.  Basic  Cathode-electrolyte-anode  construction  of  Fuel  Cell  and  electrode  reactions  [1] 


The  basic  physical  structure,  or  building  block,  of  a  fuel  cell  consists  of  an  electrolyte  layer  in 
contact  with  a  porous  anode  and  cathode  on  either  side  as  shown  in  Figure  1.  The  cells  run  on 
hydrogen,  which  reacts  with  oxygen  from  the  air  in  such  a  way  that  a  voltage  is  generated 
between  two  electrodes.  The  electrolyte  in  this  Fuel  Cell  is  the  ion  exchange  membrane  (Nation 
or  other  similar  polymers)  and  is  an  excellent  proton  (H^)  conductor.  The  electrons  flow  to  the 
external  circuit  generating  current. 

Transport  and  reaction  processes  at  the  cathode  limit  PEMFC  performance.  Oxygen  present  in 
the  air  stream  flowing  through  the  channels  in  the  so-called  bipolar  plate  (BPP),  is  transferred  to 
the  platinum  (Pt)  catalytic  particles  that  are  dispersed  in  pores  throughout  the  cathode,  through  a 
porous  diffusion  layer  and  a  porous  sub-layer  (not  shown  in  Figure  1).  Water  produced  by  the 
electrochemical  reaction  at  the  cathode  is  partially  evacuated  through  the  same  porous  medium 
to  the  air  stream. 

Limitations  of  PEM  cell  performance  at  useful  current  densities  (^0.5  A/cm^)  are  primarily 
dominated  by  mass  transfer  of  water.  The  power  density  (i.e.  W/cm^)  of  PEM  cells  is  typically 
highest  at  the  current  density  at  which  performance  starts  to  drop  drastically  due  to  reaching  the 
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so-called  “limiting  current”.  In  order  to  deploy  PEM  fuel  cells  in  high  power  density  applications 
(e.g.  automobiles  and  ships),  stable  performance  must  be  achieved  at  these  high  current 
densities.  As  the  power  density  increases,  the  weight,  volume,  and  cost  of  the  fuel  cell  system 
all  decrease,  all  of  which  are  critical  metrics  for  successful  market  penetration. 

In  order  to  obtain  high  current  density,  oxygen  must  be  transported  to  and  water  removed  from 
the  electrochemical  ly  active  catalyst  sites.  In  a  sense,  the  necessity  for  water  removal  comes 
from  its  role  in  inhibiting  mass  transfer  of  oxygen.  Since  the  PEMFC  operates  at  a  low 
temperature,  water  hang  up  (known  as  cathode  flooding)  in  a  liquid  state  is  unavoidable.  This  not 
only  impacts  PEMFC  performance  but  also  leads  to  performance  instability  [2].  Ideally,  the 
water  that  has  been  produced  at  the  cathode  is  driven  to  the  air  channel  where  the  air  would 
entrain  excessive  liquid  water.  At  the  same  time  this  airflow  would  supply  the  necessary  oxygen 
for  the  fuel  cell.  Thus,  the  directions  of  oxygen  and  water  flows  are  opposite  to  each  other. 

In  this  work  we  analyze  the  feasibility  of  trapped  water  removing  and  increase  of  the  material 
effective  permeability  by  application  of  structural  waves  to  the  system.  There  have  been  a 


Figure  2.  Structure  of  Fuel  Cell  Cathode  compartment  and  directions 
of  major  mass  flows. 


number  of  works  analyzing  waves  propagation  in  fully  saturated  porous  media  [3,  4]  -  so  called 
Biot  models.  All  known  to  us  publications  study  the  influence  of  properties  of  porous  materials 
and  fully  saturated  them  fluids  on  waves  characteristics.  In  contrast,  we  study  wave  propagation 
and  attenuation  in  partially  saturated  porous  media  and  extend  Biot’s  approach  for  multiphase 
(solid,  fluid,  and  gas)  media.  Also,  we  develop  a  novel  approach  combining  both  “elastic-wave” 
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part  and  “viscous  -  mass  transfer”  part  coupling  transport  and  acoustic  properties.  In  other  words 
we,  first,  develop  a  generic  model  allowing  us  to  evaluate  water  and  gas  transfer  induced  by 
structural  waves.  We  calculate  characteristics  of  different  waves  in  the  structure  and  couple 
them  with  water/gas  pressures  and  induced  inhomogeneous  in  space  and  time  water  saturation 
level-  S(x,t)  and,  also,  relate  it  with  oxygen  transport  answering  the  practical  question  -  could 
the  structural  waves  increase  oxygen  transport  to  catalytic  sites.  Generally,  we  develop  the  first 
model  of  acoustic  streaming  in  unsaturated  porous  structures. 


As  shown  in  Figure  2,  the  cathode  compartment  is  the  layered  structure  of  three  porous  materials 

with  different  morphology.  The  role  of  morphology  in 
wave  propagation  is  crucial  and  affects  all  model 
parameters,  such  as  Biot  parameters  and  Leverett- 
function,  in  other  words,  water  transfer  depends  on 
morphology  and  is  different  in  bi-layer  and  substrate.  At 
low  saturation  flow  is  always  channeled  and  water  can 
be  decomposed  into  connected  water  and  trapped  water. 
Trapped  water  does  not  “feel”  external  pressure  and, 
subsequently  does  not  move  to  WTP,  blocking  oxygen 
paths.  Trapped  water  is  the  “trace”  of  water  “lost  from 
the  mainstream.”  There  are  two  competitive  processes  - 
imbibitions  (from  main  stream)  and  drainage  to  the 
mainstream  and  structural  or  acoustic  waves  depending 
on  their  parameters  might  effectively  “clean”  or  vise 
versa  “spurge”  the  porous  media  changing  the  effective 
cross  section  of  oxygen  path  as  shown  in  Figure  3.  Only 
structural  vibrations  might  induce  transport  of  such 
trapped  water  “mega-droplets”.  Parameters  and 
characteristics  of  the  most  beneficial  waves  (for  example,  rotational  or  compression)  need  to  be 
defined  by  numerical  analysis  for  real-life  geometry  using  the  implemented  model. 


1^ 

Wave  propagation 

Figure  3.  Schematic  effect  of  waves  on 
water  transfer 


In  this  first  annual  report  we  present  our  results  on  executing  the  contract  according  to  the 
statement  of  work.  The  plan  of  this  technical  progress  report  is  as  follows:  1)  Problem 
Formulation  (above);  timetables  and  budget.  The  report  is  continued  with  physics-based  2) 
Model  Approach  and  Equations’  derivation  using  the  fact  that  characteristic  time  scale  of 
structural  waves/vibrations  is  much  smaller  than  characteristic  time  of  mass  transfer  and  that 
wave  amplitude  is  small.  3)  New  asymptotic  technique  has  been  developed  and  we  report  some 
obtained  analytical  solutions.  4)  In  order  to  calibrate  the  model  and  verily  the  model  predictions 
we  conduct  two  different  sets  of  experiments.  In  this  progress  report  we  review  the  first 
Experimental  Results  and  Technique  for  the  permeability  measurements  with  and  without 
applied  vibration.  We  close  this  progress  report  with  some  concluding  remarks  and  the  list  of 
tasks  that  left  to  be  executed  during  the  work  on  the  contract. 
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2. 


Scale  for  description  of  multiphase  flow  in  vibrating  porous  medium. 


For  the  description  of  the  flow  in  the  real  macro  structures  one  requires  the  scale  of  description  / 
that  is  much  less  than  the  structure  size  but  much  larger  than  the  size  of  a  pore  of  the  porous 
material  -  }.  This  allows  for  continuum  mechanics  type  of  description  of  fluid 


Liquid  F/ 


Gas  Vj 


Solid  V, 


Fig.4.  Mesoscopic  volume  V=  K5+  F/+  for  averaging  procedure.  It  is 
small  enough  to  consider  all  velocities  constant  and  large  enough  to 
include  many  pores. 


flow  and  vibration  wave  propagation,  with  the  structure  specific  being  taken  care  of  by  the 
geometry  and  boundary  conditions  of  the  domain  and  with  porous  medium  specifics  being 
incorporated  into  coefficients  of  the  resulting  system  of  equations.  The  later  must  be  either 
measured  experimentally  or  calculated  with  other  methods  that  allow  explicit  description  of  the 
material  morphology. 


The  chosen  intermediate  scale  requires  that  the  equations  to  be  formulated  through  variables 
averaged  or  smoothed  over  scale.  The  usual  procedure  (see  e.g.  [5  and  6])  is  to  split  the 

solution  domain  into  sub-domains  of  the  size  /,  Flg.4,  and  substitute  all  pertinent  quantities  like 
kinetic  and  elastic  energies,  mass  flows,  stresses  and  local  pressures,  with  the  quantities 
averaged  of  the  volume  of  the  sub-domain. 

Multiphase  materials  deformations  are  always  inhomogeneous  at  microscale,  consisting  of 


Figure  5.  Microstructure  of  porous  materials  of  FC  cathode  compartment 
various  constituents  with  different  properties.  Increasing  the  elementary  volume  of  materials,  we 
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homogenize  the  properties.  A  [minimal]  statistically  representative  material  volume  is  called 
Representative  Volume  Element  (RVE)  [6],  An  RVE  must  include  a  large  number  of  micro¬ 
elements  to  represent  local  continuum  properties.  Fuel  Cells  cathode  compartment  consists  of 
three  porous  layers  with  different  morphology,  which  predefines  the  characteristic  size  of  the 
RVE.  Thus,  for  Toray  paper  (sub-layer)  with  the  typical  pore-element  size  is  about  20  p  in  width 
and  3p  in  depth,  average  size  of  RVE  will  be  a  block  with  the  dimensions 
/x/xZ»  ~  100//xl00//x20// .  For  the  by-layer  with  the  typical  size  of  3p  the  average  RVE  size  is 

15  p  to  ensure  that  at  least  a  hundred  microstructural  objects  are  averaged  in  the  RVE  (as  can  be 
seen  from  the  micrographs  in  Fig  5.  Photo  D.  Condit,  UTRC).  In  further  consideration  all  sub¬ 
volumes  are  considered  as  points  and  all  differential  operations  are  understood  in  space  of  these 
points.  The  values  of  all  physical  quantities  are  equal  to  average  values  of  the  same  quantities 
over  the  sub-volume.  In  this  problem,  we  assign  three  characteristics  (for  solid,  liquid,  and  gas) 
for  each  material  point  instead  of  one  as  used  in  single  material  analysis.  The  following  notations 
are  used  in  each  material  point: 

|u^,u,,u^|  -  displacement  of  solid,  liquid  and  gas  fractions; 

|v^,  V,,  Vg|  -  velocities  of  solid,  liquid  and  gas  fractions; 

-  volumes  of  solid,  liquid  and  gas  phases; 


F^+V,  +  V^  =  V; 


-  volume  fractions  of  solid,  liquid  and  gas  phases; 


+S/  +  6g  1, 

-  densities  of  solid,  liquid  and  gas  fractions; 


3.  Lagrange  formulation  of  the  elastic  vibration  and  inclusion  of  dissipation 

The  complete  problem  formulation  of  elastic  waves  propagation  in  porous  medium  filled  with 
several  fluids  requires  formulation  of  conservative  part  that  results  in  vibration  of  all  constitutive 
media  and  dissipative  part  that  accounts  for  fluid  viscosity  and  elastic  energy  dissipation  in  the 
solid  matrix.  For  correct  account  of  vibration  impact  on  the  flow  dissipative  and  elastic  terms 
must  be  mixed  in  visco-elastic  constitutive  relations  of  Maxwell-Voight  type  that  involve  both 
stress-strain  and  stress-strain  rates  pairs.  Although  for  many  materials  parameters  of  visco¬ 
elastics  constitutive  relations  have  been  measured  and/or  computed,  the  same  for  porous 
materials  and  especially  for  the  recently  invented  porous  membranes  are  absent.  Here,  we 
suggest  simplified  approach  that  is  based  on  the  fact  that  the  time  scales  of  vibration  induced 
with  the  external  source  in  kHz  range  and  higher  and  filtration  of  fluids  in  the  porous  matrixes 
are  very  well  separated.  Additionally,  computational  implementation  becomes  much  more 
simple  as  it  allows  natural  computation  of  elastic  part  in  Lagrangian  coordinates  and  filtration 
part  in  the  Eulerian  ones. 
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The  problem,  therefore,  is  to  split  the  formulation  of  Free  Energy  functional  for  the  whole 
problem  into  Lagrange  function  for  elastic  vibration  of  all  constituents,  which  has  originally 
been  formulated  by  Biot  [3],  and  the  part  that  is  responsible  for  filtration  due  to  slow  varying 
gradients  of  hydrostatic  pressures.  Subsequently,  various  dissipation  mechanisms  are  included  to 
obtain  generalization  of  the  Darcy  law  and  hysteresis  phenomena. 

Thus,  the  goal  is  to  derive  expressions  for  kinetic  energy  K,  elastic  energy  t/and  dissipative 
function  R  in  terms  of  averaged  over  representative  volume  physical  quantities.  Having  in  mind 
that  generally  all  three  functions  depend  on  {u,  v.au/ac},  the  final  equations  have  the  following 
form: 


d 

^  dL  ^ 

d 

_ 

(  dR  ] 

1 

_^_dR 

dx 

^d{du/dx)  ^ 

dx 

^d{dvldx)^ 

5/1 

^dvj 

5u  dv 

where  L  =  K-U  and  x  is  coordinate  vector. 


Theoiy  of  liquid  flow  through  porous  media  uses  equation  of  mass  balance  closed  with  the  help 
of  Darcy  equation.  It  is  the  purpose  of  this  work  to  show  how  vibration  waves  described  by  the 
system  (1)  influence  filtration  rate.  The  main  philosophy  of  this  approach  is  to  add  small  elastic 
vibrations  and  mass  transport  induced  by  them  “on  top”  of  viscous  flow.  It  is  done  in  the  further 
paragraphs. 


3.1  Kinetic  energy  K. 

Original  Biot  equations  [3]  have  been  derived  for  only  one  phase  filling  the  porous  medium. 
Two  or  more  phases  coexisting  inside  porous  matrix  require  additional  treatment  and  general 
pidelines  on  how  to  write  correctly  the  mass  tensor.  The  correct  expressions  for  the  mass  tensor 
in  fully  saturated  by  single-phase  porous  media  (Biot  problem)  have  been  recently  obtained  by 
an  averaging  procedure  in  [7].  We  generalize  this  approach  for  the  multiphase  situation  in  order 
to  obtain  consistent  expression  for  the  kinetic  energy  and  inertia  force  as  functions  of 
constituents. 

We  are  looking  for  expression  for  kinetic  energy  K  =  +  p,\]  +  /j^v^  )  through  averaged 

over  some  small  volume  of  porous  medium  V  values  of  ((v,),(v,),(vg)),  where  symbol 
{•••)  =  ^  |^^(•••)  "leans  averaging  procedure  over  the  RVE  as  shown  on  illustration  Figure  5. 

Each  phase  velocity  |vj,v,,Vg }  can  be  re-  written  as  a  sum  of  its  average  value  and  a  deviation 
as  follows:  |(v,)  + Av,,(v,)  + Av,,^v^^  + Av^}.  Spatially  averaged  kinetic  energy  is  a  sum  of 
terms  with  averaged  squares  of  partial  velocities,  which  can  be  expressed  in  the  following  form: 
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W  =  {(A^')  +  A(v/)  +  Pg(v5))  = 

^Ps^s  (({v. )  +  ^  +  p,s,  ^(( V, )  +  Av;  )^  ^  +  PgSg  ^(( Vg  )  +  AVg  )  ^  j 


(2) 


It  is  important  to  note  that  volume  fractions  appears  before  averaging  operation  due 

to  the  fact  that  velocities  |v^,v,,v^|  are  non-zero  only  in  their  respective  sub-domains 
I  within  a  RVE.  Modifying  (2)  we  get  the  expression  for  the  averaged  over  RVE 


kinetic  energy: 

{K)  (v,  f  +  Pl^l  (v;  >^  +  Pg£g  (vg  }^  +  Ps^s  f  )  +  Pl^l  ((^V/  f  "j  +  PgSg 


Due  to  solid  matrix  rigidity  the  term  Ps£s{{^"^s)l  usually  small  and  we  neglect  it  from 

further  consideration.  To  close  system  of  equations  of  motion  one  has  to  assume  relations 
between  |(J(Av,)^^,^(AVg)^^}  and  {(v,),(v,),(vg)}  as  was  pointed  out  in  [7]. 

To  obtain  classical  Biot  kinetic  energy  terms  for  two  phases  flow,  only  ||^(Aw ,  ^(aw  ^  ^  | 

should  be  expressed  through  mean  velocities  of  constituents: 


^  figi 


This  parameterization  is  similar  to  the  assumption  of  Poissonian  statistics  for  velocities 
fluctuations.  The  final  expression  for  the  average  kinetic  energy  in  this  case  is  as  follows: 


2{K)  =  p,s,{y,f  +p,€,{y,f  +p^s^{y^f  +{p,€,fi,  +  Pg£gfig,)({y,)-(ys)f  + 

[Pi^iPgi  +  Pg^gPg)[{^g  )  -  {^s  >) 


One  phase  situation  that  has  been  treated  by  Biot  corresponds  to  only  one  term  -^(Av,)^^  and 
2p,+\  =  a-  tortuosity  in  Biot’s  notation.  In  the  two-fluid  model  describes  correlation 
between  relative  motion  of  gas  versus  solid  and  correspondent  motion  liquid  and  vise  versa. 


The  dynamic  force  acting  on  a  system  is  the  derivative  of  the  kinetic  energy  with  respect  to 
velocity.  Using  (4)  we  can  obtain  the  following  form  for  (inertia)  D’Alamber  force: 
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It  is  important  that  all  elements  of  mass  matrix  explicitly  depend  on’  volume  fractions  of 
constituents,  that  easy  to  relate  with  local  fluid/gas  saturations.  This  is  the  new  result  because 
classical  Biot  theory  considers  only  fully  saturated  porous  media  and  the  solid  body  porosity  is 
the  primary  characteristics  of  that.  Now,  the  mass  matrix  depends  on  volume  fractions  and  cross- 
correlations  between  relative  phases  motions  that  need  to  be  approximated  and/or  experimentally 
measured.  Usap  of  the  averaging  procedure  allows  explaining  of  non-diagonal  form  of  the  mass 
tensor  in  the  Biot  model  [3].  It  also  provides  connection  between  experimentally  measured  fluid 
densities  and  effective  densities  that  is  to  be  used  in  equation  for  wave  propagation  (I). 

3.2  Elastic  energy  U 


Usage  of  averaging  procedure  of  the  previous  paragraph  allows  explaining  non-diagonal  form  of 
the  mass  tensor  in  the  Biot  model  [3].  It  also  provides  connection  between  experimentally 
measured  fluid  densities  and  effective  densities  that  is  to  be  used  in  equation  for  wave 
propagation  (1).  One  can  also  try  to  use  the  same  averaging  procedure  as  was  used  above,  for  the 
derivation  of  the  averaged  elastic  energy  U  as  a  function  of  displacements  averaged  over  RVE. 
This  would  be  analogous  to  derivation  of  turbulent  models,  fluctuations  being  imposed  by  the 
randomness  of  the  porous  matrix  rather  than  by  the  local  flow  instability.  However,  as  numerous 
attempts  to  derive  elastic  constants  show  [4],  those  expressions  strongly  depend  on  material 
morphology.  Here,  we  follow  more  pragmatic  approach  of  Biot  [3]  and  derive  expressions'  for 
elastic  energy  with  parameters  being  measurable  in  specified  set  of  experiments  or  calculated 
using  direct  microscale  modeling.  More  important  that  for  our  purpose  is  sufficient  to  combine 
correctly  equally  phenomenological  theories  of  liquid  stagnation  in  porous  media  with  the  Biot- 
type  vibration  equations. 

With  standard  notations  of  theory  of  elasticity  small  strain  is 


duj  j 

+  "7  ^  =  ^11+^22+  ^33  *  is  dilatation.  Because  there  are  no  stresses  in  the 

ax,,  j 

elastic  body  at  zero  strain,  and  stress  is  the  derivative  of  elastic  energy  with  respect  to  strain, 
there  should  be  no  linear  terms  in  the  energy  expansion  on  strain.  Terms  with  power  higher  then 
three  would  mean  existence  of  several  metastable  states  typical  for  phase  transformations  and 


‘  Of  course,  we  use  phenomenological  considerations  here,  so  it  is  derivation  in  physical  but  pure  mathematical 
sense. 
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not  existing  in  elasticity.  The  expression  for  elastic  energy  can  be  written  in  the  general  form  as 
follows: 


-X-6^  +2u-'y  e  ■  e„  +  self  elastic  energy  of  porous  solid  (5) 

ry  S  *  IJ  Ij 

^  y 

3/  {9^,6i)  +  ^g[9^,0g^  +  energy  of  liquid  and  gas  dependent  on  solid  distortion 

We  assume  there  is  no  shear  in  gas  phase.  Additionally,  taking  into  account  that  liquid  flow  rate 
in  porous  media  is  very  low;  we  may  neglect  the  shear  component  in  liquid  strain  tensor. 
Therefore,  the  elastic  energy  for  gas  and  liquid  phases  depends  only  on  volume  change  or,  in 
other  words,  strain  dilatation.  Stress  tensor  for  solid  matrix  and  pressures  for  liquid  and  gas  can 
be  obtained  by  taken  correspondent  derivatives: 

»  se,  *  se. 

„  s  ac/  (6) 

.  du 

«  gg^  gg^ 

Here  <5'^  is  Kronecker  delta  with  components  of  unity  matrix.  Functions  account  for  the 

pressure  change  due  to  liquid  and  gas  motion  under  solid  deformation.  The  models  for  3,  ^  {O^ ,  O,  ^  )  are 
discussed  in  separate  paragraph.  Below  we  consider  some  special  cases  and  the  corresponding  potential 
energy  forms. 


3.2.1  Small  vibrations 

In  the  absence  of  large  motions  and  near  equilibrium  between  all  three  phases,  U  can  be 
expanded  up  to  the  second  terms  in  ^^,9, g).  As  has  been  already  discussed  above,  it  means  that 

at  this  internal  state,  reaction  of  all  phases  on  applied  deformation  is  reversible  and  elastic.  We 
emphasize  here  that  the  theory  we  are  constructing  in  this  work  is  elastic  correction  to  the  “pre¬ 
existing”  viscous  flow.  This  assumption  is  fully  justified  for  small  vibrations,  when  signal 
magnitude  does  not  exceed  5%  of  structural  characteristic  size.  Naturally,  bigger  loads  would 
lead  to  the  structure  fatigue  failure. 


Thus,  in  the  absence  of  large  motions  and  near  equilibrium  between  all  three  phases,  (/can  be 
expanded  up  to  the  second  terms  in  {^,,0,^ ),  where  (G,^  )  is  defined  as 


where 

-e-j  -  self-elastic  energy  of  porous  solid; 

^  0 

^si  =Mis -0^  •6>-  energy  of  interaction  between  liquid  and  solid  and  gas  and  solid; 

V,  -0]  -self-elastic  energy  of  liquid  and  gas; 

‘01  ‘0g  -energy  of  interaction  between  liquid  and  gas; 

-0,  +  -0^  -energy  of  interaction  between  liquid  and  gas  and  correspondent 

static  pressures;  and  simple  form  for  3/  ^  =  P,,^  -s,^  has  been  chosen. 


Here,  M,  and  are  self  influence  bulk  coefficients  for  liquid  and  gas  correspondently, 

^i>  O' ~  g/)  are  coefficients  of  mutual  interaction.  Constitutive  relations  relating  stresses 

(pressures)  and  strains  are  obtained  by  differentiation  of  elastic  energy  with  respect  to  strain 


tensor  and  take  the  form: 


n>,  =  -s,-p, -S,  =  — =  M,  •«,  +  M„  e,  +A/^ 


3.2.2  Biot  theory 

In  his  original  paper  [3]  Biot  considered  fully  saturated  porous  media  with  only  one  liquid.  He 
also  disregards  porous  pressure  P,'’"" .  Under  these  simplifications,  from  relationship  (7)  setting 


(Here,  ^  0^,  0^^  are  quasi  equilibrium  values  of  liquid  and  gas  dilations  caused  by  static 

pressures,  conditions  of  quasi  equilibrium  being — ^  =  Q .) 
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all  gas-related  parameters  as  well  as  zero  we  obtain  the  constitutive  relations: 

cr,5  =  =  A  •  +  2  •  •  gy  +  M,  •  7;  •  6>, 

(8) 

<y\j=-<l>-p,-Sy=^  =  Mre,+MrYres 

We  substituted  liquid  volume  fraction  from  (7)  with  material  porosity,  which,  by  definition, 
equals  volume  fraction  of  pores  in  the  solid.  Due  to  the  fact  that  solid  is  fully  saturated;  liquid 
volume  fraction  is  identically  equal  to  the  porosity.  Thus,  we  deduced  Biot’s  constitutive  model 
from  ours  as  a  limiting  case. 

By  this  point  we  have  considered  Lagrangian  and  its  constituents.  In  other  words,  we  have 
derived  all  relations  describing  the  elastic  behavior  of  the  system.  In  order  to  complete  “elastic 
part”  modeling  and  before  we  switch  to  dissipative  effects,  we  would  consider  developing 
methodology  for  constitutive  model  coefficients  ( M,  and  Yi  )  determination. 


4  Elastic-Biot-type  coefficients 

4.1  Effective  Media 


First,  all  Biot  coefficients  can  be  expressed  directly  from  the  effective  bulk  moduli  of  porous 
media,  matrix,  and  a  fluid  [4].  The  problem  of  estimating  the  effective  properties  is  not  trivial  but 
might  be  performed  by  methods  developed  in  the  theory  of  mixtures  [8].  Generally,  the  best  we 
can  do  is  to  predict  upper  and  low  bounds  of  the  effective  parameters,  because  its  precise  value 
depends  on  geometrical  details.  It  is  important  to  note  that  all  such  estimates  work  well  only  for 
isotropic  media.  The  extended  review  of  methods  is  given  in  [8]. 


The  best  achievable  bounds  without  specifying  geometrical  details  are  the  Hashin-Shtrikman 
bounds  [9]; 

j^HS±  _ _l _ ^2 _ . 


(  +  2/2|  ) 


where  K;  are  bulk  moduli  of  phases,  pi  are  shear  moduli  of  phases,  and  e;  are  volume  fractions. 
These  estimations  for  Toray  paper  (in-plane)  are: 


1 .2  •  10^  Gpa  =  17.4  Psi  =  <K<  =  3.33  Msi  =  23  GPa ,  where  initial  fibers 

compressibility  has  been  taken  equal  20 Msi— \30GPci  Mjjher  —25msi=l60GPa 

Bounds  for  shear  modulus  are  0  s  <p<  =  4.5  Msi  =  31  GPa.  Difference  in  upper  and 

low  bounds  are  extremely  large  as  expected,  because  there  are  significant  differences  in  matrix 
and  inclusion  properties.  Experimental  data  for  in-plane  elastic  constants  varies  from  40  Ksi  till 
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140  Ksi.  We  assume  that  K,^^,  =60  Ksi^.  We  apply  empirical  approach  of  Marion  [8,  p.  177] 
and  calculate  the  weight  coefficient,  which  reduces  the  upper  bounds  to  experimentally 

measured  values:  w  =  « 0.018.  Thus  we  should  expect  in 

M  —M 


plane »0.018-31  GPa- 558/fPa  =  80  AIs/.  Measurements  throughout  the  thickness  gives 
extremely  low  number  of  «1500;?5/.  This  might  be  modeled  using  Hashin-Shtrikman 
bounds  if  we  imagine  that  fibers  contact  only  over  small  “point-like”  area  and  effective  porosity 

in  this  direction  would  be  «  0.9997 .  Indeed,  £• « 1  -  =  j  =  0.9997  and 

/  ^between 


«5000  psi\ 

There  is  also  known  test  result  that  compression  modulus  of  the  composite  of  Toray  paper  and 
bi-layer  together  is  1699  psi.  Assuming  weights  proportional  to  the  thickness  of  each  layer  (25 
pm  and  175  pm)  we  might  estimate  « 1.5  -2.5  Ksi.  Hashin-Shtikman  estimates  based  on 


dual  porosity  give  bounds  150p^/  =  ^K<  =  7Msi  =  50GPa.  Such  a  big  difference 

arose  because  we  did  not  take  into  account  the  compactions  of  sintered  powder.  As  noted  by  M. 
Kachanov,  the  stiffness  of  the  media  decreased  when  the  pore  shape  is  concave  and  we’d  better 
use  Mori-Tanaka  approximation. 


4.2  Estimation  of  Biot  Coefficients  for  Solid-Liquid-Gas  System 

The  “elastic”  system  that  we  consider  now  has  the  following  general  form: 

^y=-^i'Pr^ij~^h‘^s'^^r^i+^gr^g  (10) 

First,  the  dependence  of  the  Biot  coefficients  on  saturation  is  evaluated.  It  is  clear  that  to  satisfy 
the  relations  above  near  the  limits  £/  ->0  and  >  0,  the  elastic  constants  must  have  the 
following  dependency  on  saturations: 

^Is  =  ^r^is'-’^gi  =  £i  ’mg,;M,  =  e, -nii 

•  f»gs'Mg,  =^g-^gi-Mg=£g-mg 


From  this  it  follows 


’  We  understand  that  there  were  attempts  to  measure  Young  modulus,  so  Bulk  modulus  should  be  about  3  times 
smaller  for  isotropic  materials. 

Porosity  here  depends  on  the  area  of  contacts  and  might  vary  a  lot.  We  also  can  approximate  the  experiments  if  we 
use  previously  described  Marion  weight  factors  and  d  jm ,  and  subsequently, 

^  "  %00  *  ^  0  01 8  •  =  75  Ksi)  =  1 350  psi. 
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Mg,=£r^g-mg, 

Comparing  Biot  coefficients  [3]  for  two-phase  case: 


where  index  /  stands  for  “fluid”  (gas  or  liquid),  with  the  equations  for 
[M,g,M,sMgsMg,}  above  one  gets: 


(12) 


And  finally,  defining  moduli  ^  ;  Kj  =  one  gets: 

{K,-K,y0,+M,,-O,+M^,-0g=O 

M,s-0s  +(^1  -£iPi)-0i  +Mgi -^g  =0 

Mg,-0,  +M^, -0,  +(Mg -SgPg)-0g  =  0 

This  homogeneous  linear  system  has  a  non-trivial  solution  only  if  its  determinant  equals  to  zero: 

K,-K,  M,  M,, 

M„  M,-£iP,  M^,  =0  (15) 

^gs  ^gi  ^g-^gPg 

Solving  this  equation  with  plausible  assumption  that  Kg  «  Kj  «  ,  and  substituting 

expressions  for  coefficients  above,  we  obtain 
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4.3  Simplified  Expressions  for  System  Coefficients 

Consider  stiff  matrix  and  assume  that  the  equalities  «  K,  «  and  K  «  are  satisfied. 
In  this  case,  coefficients  have  very  simple  form: 

■Af/=f/Ar,;  ^gs~ - --Eg-Kg,  A/^  =£■  A”  ;  and  finally. 


a:.  k  ]\  £ 


•S,-Sg. 


More  accurate  expression  for  the  interaction  coefficient  is 


Mg,. -2. 


ErSg 


It  can  be  simplified  to  the  form  given  above  or  neglecting 


the  term  ,  to  the  simple  expression  M„,  «  -2—^ 
K,  K 


KiKg(\- 


■e,-€g. 


Thus,  for  Toray  paper  the  full  set  of  constants  is  =  70  ksv,  =  35  ksi\  =  0.75;  for 

through  thickness  we  need  to  use  =1.2  ksv,  =  0.7  ksi  and  the  value  of  porosity  in  this  direction 

is  not  clear  at  this  time.  Density  of  Toray  paper  is  p^  =450-1 634  .  For  the  bi-layer  data  are  very 

few  and  far  between.  However,  the  first  shot  might  be  =  1  -  2  ksv,  //^  =  1 .5  -  3  ksv,  (j)  «  0.25;  and 

density  =  1 634-^ .  For  further  use  one  needs  the  following  data: 

^ water  -  2-2  GPa;  =  1 .2  •  10"^  GPa .  Viscosity  of  water  is  =  1 0”^  Pa  ■  sec  and 
=10“^ Pa  sec. 

Closing  this  section  we  emphasize  that  here  we  report  novel  results  revealing  dependence  of  the 
constitutive  model  coefficients  from  the  species  volume  fractions  or  saturation.  We  recognized 
and  derived  exact  expressions  showing  that  dynamic  stress  -  strain  relationship  strongly  depend 
on  local  saturation  of  porous  media.  Thus,  the  pressure  distribution  under  acoustic  signal  is 
naturally  inhomogeneous  in  an  inhomogeneously  saturated  media. 


We  have  completed  the  analysis  of  the  “elastic”  deformation.  Next,  we  focus  on  “viscous- 
dissipation”  part,  and,  then;  couple  them  in  a  unified  model.  Inelastic  effects  are  added  to  the 
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model  by  addition  of  appropriate  dissipative  function.  This  function  will  describe  flow 
phenomena  and  solid  stress  attenuation. 


5.  Dissipative  function  R 


Dissipative  function  determines  amount  of  energy  that  mechanical  system  losses  per  unit  time 
[10].  Dissipative  processes  are  always  irreversible,  so  the  dissipative  function  should  not  change 
sign  and,  consequently,  to  the  first  approximation  should  be  positive  definite  quadratic  form.  We 
assume  there  are  three  major  dissipative  processes  in  the  porous  partially  saturated  continuum. 
First,  it  is  energy  losses  associated  with  the  friction  due  to  liquid/gas  phase  relative  flow  with 
respect  to  the  solid  matrix.  We  would  call  this  D’Arcy  term  ,  emphasizing  the  mass- 

transfer  role  in  this  process.  The  second,  energy  is  lost  due  to  waves  damping/attenuation  in  the 
solid  matrix,  or  in  other  words,  terms  responsible  for  viscous  friction  in  dry  solid  matrix. 
The  third  term  is  the  energy  dissipation  rate  due  to  spatial  gradients  of  relative  solid-liquid 

and  solid-gas  motion.  This  term,  in  particularly,  describes  hysteretic  effects. 


R  =  R^'^^+R^‘‘‘+^.= 


f  K^'  (w, ) + f  K-'  (wg ) 

^  •  (-9.  )^  +  2  •  77  •  (ta/ ,  <Wg ,  5, ) 


Here  we  denote  throughout  vectors  w,  ^  =  £,  g  ■  (v,^  -  )  filtration  velocities  of  liquid  and  gas. 

We  also  introduce  symmetric  part  of  velocity  gradient  tensor  for  solid  matrix  -  strain  rate  tensor 
and  solid  strain  dilatation  rate  -  i9^  as  well  as  filtration  dilatation  rate  for  liquid  and  gas  as  follows: 

1  fsv''  1  ,  5  c  dw‘{^  dw2^ 

\  j  '  j 

In  fluid  mechanics  equations  velocity  ( v  =  m  )  plays  the  same  role  as  displacement  (  m  )  in  solid 
mechanics. 


Dissipation  function  describes  inner  friction  and  should  be  zero  if  there  is  no  internal 
motion  [10],  particularly  in  case  of  rigid  motion.  It  means  that  this  term  should  depend  on 
velocities  gradients.  For  isotropic  body  two  viscosities  (bulk  and  shear)  should  be  introduced  as 
it  shown  in  (18). 

In  special  case  of  weak  liquid-gas  interaction  ‘S{g,{cOi,(Og,3^)  can  be  simplified  to  just  a  sum 
Functions 91/^ (<y,g)  and  are  discussed  in  special  paragraph 

below.  Under  assumption  of  small  oscillations  near  equilibrium,  [a, ,0)^,3^)  has  a  general 
quadratic  form  [11]: 
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^A^n<»gA)=^R, -(of  -al  +  R^i-o), -eo^+R,,  o),  •  .9,  + /?^^  •  •  ,9, 

In  this  case,  the  corresponding  friction  forces  are  linear  with  filtration  dilatations  rates  and  equal 
to  zero  at  zero  dilatation  rates.  To  include  hysteresis  into  this  formalism,  non-quadratic  form  of 
*^ust  be  preserved. 

According  to  the  basic  equation  (1)  the  corresponding  dissipative  forces  can  be  obtained  from 
the  dissipative  function  by  derivation  with  respect  to  appropriate  velocities  and  velocities 
gradients.  We  start  with  forces  caused  by  filtration.  The  first  three  terms  in  the  right  hand  site  of 
(18)  give  the  general  expression  for  Darcy  type  forces: 


^n'Arcy  _  dR 

dv^  dv,  dVg 

dy,  K. 


Mr  Mg 


*  dv^  ^  K,  ^  ^ 


^i^g. 


As  eveiywhere  in  this  report,  subscripts  (s,  I,  and  g)  correspond  to  solid,  liquid,  and  gas  phases 
accordingly.  K, ,  K^,  and  K^^are  tensors  for  liquid  permeability,  gas  permeability  and  cross  liquid- 
gas  permeability.  Last  term  characterizes  the  flux  of  one  of  species  when  pressure  is  applied  to  another. 

All  permeability  tensors  components  are  functions  of  liquid  saturation  5  =  orf,  and  accordingly. 

The  velocities  gradients  dependent  terms  contribute  to  internal  stresses  and  need  to  be  added  to 
expressions  (7)  for  completion  of  constitutive  relations: 


- - g-&,-6,j+2-r]-v,j-e, 


/  diss  _  diss  o  (^/ )  « 

-srP,  =  (20) 

diss  _  .  p  diss  ^  „  ^g  {^g  )  y 

For  filtration  part,  an  assumption  is  that  //,  are  the  same  as  those  for  acoustic  part.  In  other  words, 

frequency  dependence  of  is  neglected.  Thus,  relaxation  forces  for  filtration  of  gas  and  liquid 
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^■"X'"^‘^=-^  =  {^,)  ^K7‘V,+(^,) 


k;;v. 


filtrjpD'Arcy 


Here,  absolute  velocities  of  liquid  and  gas  V,  ^  stand  instead  of  relative  velocities 

w,  g  g  •  ( V/  g  -  )  as  filtration  velocity  of  solid  is  zero  -  non-consolidating  porous  material. 


saturation  is  defined  as 


6  J-Leverett  function,  dissipative  functions  and  hysteresis. 

To  provide  a  model  of  two  non-mixing  fluids  that  includes  effect  of  residual  saturation  we  have 

to  account  for  strong  dependence  of  material 
Llq»l<l.air  Interfaces  in  I  properties  on  saturation  and  wettability.  Liquid 

^kine  of  rods  (Collins  saturation  is  defined  as  S  =  ^=-^. 

it  _ I  ^  1-f^ 

Capillary  pressure  Pc  is  usually  defined  as  the 
T  \  ^  \  pressure  difference  between  two  (wetting  and 

"  \  ’  non-wetting)  phases  as  a  function  of  the 

\  „  „  saturation 

5'*=  ~  ° 

^“‘^0  Therefore,  we  define  the  capillary  pressure  as: 

- —  p  (S)  =  p.-P  For  modeling  and  correlation 


100%  80 


iJlc 

Figure  6.  Schematic  explanation  of  J-Leverett  function  .  ^ 

lary  pressure  can  be  described  by  a  dimensionless  | 

,  .  ,  APiS)  [k  ^  •  “  4 

Leverett-J  function:  J  =  — —J—,  where  a  is  |  ^ 


interfacial  tension  between  the  two  phases  and  ^the 
porosity  of  the  solid  matrix.  Definitely,  we  can  write 
down  the  same  kind  of  functions  for  each  phase  and 
their  difference  would  determine  the  capillary 
pressure.  In  our  modeling  we  would  call  J,g{£,g) 

Leverett  function  for  liquid  and  gas  separately.  In 
principle  they  depend  on  as  well,  but  the 

dependency  is  unknown  and  requires  additional 
measurement  and/or  microscale  simulations.  One 
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Figure  7.  Hysteresis  in  drainage-imbibitions 
process  (From  Bear  [15]) 
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can  identify  with  Leverett  function  correlates  data  over  several  orders  of 

magnitude  change  in  permeability  as  stated  in  [15].  Capillary  pressure  is  a  measure  of  the 
curvature  of  the  interface  separating  phases,  therefore  Leverett  J-function  is  a  strong  function  of 
morphology.  During  a  subsec]uent  water-drive  and  turns  off  some  amount  of  water  eventually 
remains  capillaiy  trapped,  floating  as  disconnected  blobs  in  the  pores.  The  residual  water 
saturation  So  is  determined  by  the  pore  space  morphology  and  the  structure  wettability.  Sq  is 
called  the  “irreducible  saturation”,  which  correlates  with  the  amount  of  trapped  water  and  for 
wetting  phase  appears  to  be  around  0.08  [16]. 


It  is  important  to  note  that  J-Leverett  function  describes  quasi  -  equilibrium  state  and  small 
deviations  from  the  equilibrium,  which  corresponds  to  slow  species  flow  through  the  porous 
media.  As  has  been  already  mentioned,  this  slow  flow  does  not  appear  in  Biot-type  dynamic 
(wave)  equations.  Only  combined  “slow”  and  “fast”  motions  describe  mass  transfer 
phenomenon.  General  expression  of  the  appropriate  (Leverett)  term  of  the  Lagrangian  or  free 
energy  is  the  function  of  the  saturation  and  material  properties^  Dropping  all  details,  for  the  sake 
of  clarity,  we  can  say  that  in  the  first  approximation  (S)  +  (S) .  This  form 


leads  to  the  expression  for  capillary  pressure  P^(S)  =  P,-P  = 


QRfast 

-  +  -^  =  J{s)  +  F{S). 

Slow  fast 


Dependency  of  constitutive  relations  on  5  =  -^  leads  to  the  hysteresis  [12].  The  rate  dependent 


part  might  be  specified  in  Lagrangian  by  9?^  ^  ^  ^  *  Fj  ^  {^)dg .  In  general,  it  is  multi¬ 

valued  function  with  even  power  of  terms  to  insure  it  positive  definiteness.  Piecewise-linear 
approximations  have  been  suggested  [12]. 


As  shov^n  on  illustration  in  Figure  7  for  wettable  materials  to  drain  fluid  from  the  porous 
material  the  higher  pressure  should  be  applied  then  to  imbibe  to  the  same  saturation.  It  means 
that  when  we  have  pressure  oscillation  in  a  hydrophilic  material  the  residual  water  should  stay  in 
the  porous  media  and  vise  versa  for  hydrophobic  materials.  The  first  conclusion  from  this 
qualitative  picture  we  can  do  is  that  one  need  to  use  hydrophobic  substrates  in  fuel  cells  to 
appreciate  cleaning  effect  of  pressure  oscillations. 

7.  Equations  of  filtration  for  two  non-mixing  fluids  and  generalized  D’Arcy 
law. 


Fluid  filtration  in  porous  media  is  traditionally  described  by  equations  based  on  mass 
conservation  [5, 12,  13,  and  15]  rather  than  on  Biot  type  equation  [3]: 

!■(«,/>, )+V(p,(V,+w,))  =  0 


Recent  advances  in  thermodynamics  of  multiphase  flow  in  porous  media  [for  example,  5,11,12]  provide  a 
theoretical  framework  for  a  proper  and  consistent  definition  of  Leverett  function  in  acoustic  field. 
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In  the  case  of  saturated  incompressible  fluids  those  equations  are  simply  direct  consequences  of 
the  definition  of  the  dilatation;  =  V  •  (u^  -  u, ) .  In  the  two-fluid  case  when  saturations  and 

change  both  in  time  and  space,  situation  becomes  more  complicated  especially  taking  in  mind 
that  Leverett  functions  and  are  measured  as  functions  of  £■/  and  rather  than 

9,g-^,g\p,,£,g)-  One  can  consider,  therefore,  the  two  continuity  equations  as  equations  for 
additional  variable  £,g  and  solve  them  together  with  Biot  equations  for  vibration.  Note  that 

because  the  continuity  equations  are  conservation  laws  and  expressed  as  full  time  derivatives, 
their  contribution  to  Lagrange  equations  taken  with  the  Lagrange  multipliers  don  t  change 
equations  of  motion. 


Note  that  when  liquid  and  gas  can  be  considered  incompressible,  densities  can  be  cancelled  in 
(21)  and  rate  hysteresis  function  from  the  Leverett  expression,  which  in  general  depends  on  flux 
dilatation,  can  be  simplified  to  the  form  mentioned  in  the  previous  section  as  follows: 


r 


hg 


r 


dt 


=  F, 


i,g 


dt 


(22) 


Equilibrium  condition  requires  that  the  resultant  of  all  applied  forces  be  compensated  by  internal 
stresses.  For  the  liquid  the  total  pressure  gradient  should  be  equal  to  the  sum  of  D’ Alamber  and 
D’Arcy  forces  (see  equations  7  and  10). 


V  (-f,  Pi )  =  (23) 

Next,  without  loss  of  generality,  we  neglect  all  “gas  terms”  and  obtain  the  simple  form  of  the 
generalized  D”Arcy  law  for  unsaturated  porous  media  under  acoustic  loading.^  Expressions  for 
liquid  components  of  these  forces  are  given  in  the  derived  expressions  (4a)  and  (19)  accordingly. 
Note,  that,  as  everywhere  in  this  work,  hydrostatic  pressure  has  the  sign  minus  in  accord  with 
continuum  mechanics  sign  agreement.  Thus,  equation  (23)  can  be  re-written  in  the  following 
form; 

V {-£iPi )  =  ( V/  -  V, )) + ^{PiF^i ) = 

The  D’ Alamber  term  depends  on  both  relative  (vi-v*)  and  absolute  liquid  velocities.  Resolving 
this  equality  with  respect  to  the  flux  q  =  w ,  and  also  taking  into  account  that 


Pi  =  -M6i  -My6^,v^q  immediately  obtain: 


Pi 


dt 


5(u, -uj 


dt 


Pi  5 
+  — T-l 

Si  dt 


5u, 

dt 


(24) 


Here  the  flux  tensor  depends  on  all  components  gradient  of  pressure,  dilatations  and  kinematical 
terms.  This  is  the  novel  generalization  of  the  D’Arcy  Law  for  the  acoustically  loaded 
unsaturated  porous  media.  The  first  term  represents  generalization  of  traditional  pressure 
gradient  corrected  for  the  variable  saturation.  Very  important,  that  saturation  gradient 


*  In  all  further  derivations  and  calculations  we  use  the  full  system  and  simplified  expression  here  only  for  illustrative 
purposes. 
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P  ^  -  PVltif  is  also  strong  driving  force  specific  only  for  unsaturated  inhomogeneously  filled 

porous  structure.  This  term  plays  an  important  role  but  only  if  either  external  pressure  or 
external”  dilatation  has  been  applied.  The  second  term  explains  the  role  of  solid  and  liquid 
deformation  as  a  driving  force  for  the  flux.  Due  to  dilatation,  fluid  is  squeezed  out  of  the 
compressed  pores  due  to  induced  increased  fluid  pressure  in  the  “extended”  area  with  reduced 
liquid  pressure.  The  two  last  terms  show  the  role  of  dynamics  (acoustics)  on  the  liquid  flux. 


Using  the  new  formulation  (24)  and  assuming  everywhere  that  saturation  does  not  change 
significantly,  we  might  estimate  orders  of  magnitude  for  the  required  parameters  of  structural 
acoustics  impact  to  achieve  the  comparable  effect  with  static  pressure.  D’Arcy  term  for  the  fuel 
cell  cathode  compartment  is  generated  by  application  of  the  pressure  of  1  psi  over  the  distance  of 

'P^>  =34£6^^.Utting 


200  m.  Thus,  the  “static”  pressure  gradient  is  of  the  order  of  Vp : 


200  pm  m 

magnitude  be  1%  of  the  thickness,  we  immediately  obtain  that  a  dynamic  term  has  the  order  of 

•  A  2  2 

magnitude  p^  —  ^pAa  »l.£-3<y  and,  subsequently,  equating  these  two  pressures,  we 
obtain  that  the  required  frequency  is  about  30  KHz. 


Constitutive  reiations  and  compiete  system  of  governing  equations 


By  this  point  we  have  described  all  forces  acting  in  the  system  and  the  system  possible  reactions. 
We  introduced  elastic  and  dissipative  models.  Combination  of  “elastic”  and  “viscous”  equations 
completes  physics-based  system  and  allows  us  to  predict  the  influence  of  vibration  on  mass 
transfer  and  dynamics  of  saturation  on  vibration  parameters.  Using  these  derived  partial  models 
one  can  get  the  following  constitutive  relations: 


,  scfiss 

(Ty+ay 


eu  dR 

dey  dVy 


=  ^-0s-Sy+2-p-ey+g&^-Sy+2-r]-Vy+M,,-e,+M^,-e^ 

-^rFiMSy-£gF^[<a^)Sy 


(Ty+CTy 


du  ^  m 

dOf  dco, 


=  •  [p,  +  \Sy=M,-e,+ +  M^, 


■Og  +£, 


■FiM 


(25) 


gdiss  _  dU  d9{ 


d0„  do)„ 


— ^*0 


K) 


Here,  each  of  the  matrix  equations  corresponds  to  solid,  liquid,  and  gas  phases  accordingly.  As 
was  mentioned  already  7  is  a  static  part  of  Leverett  function  for  each  component  (liquid  and  gas) 
and  is  a  hysteresis  part.  These  constitutive  relations  should  be  supplemented  by  the 
definition  of  the  Leverett  function,  simplified  form  of  which  is: 
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(26) 


ns,y^-p,+p,=j<.s,) 

Some  particular  expressions  for  7(s)  and  J(£)  can  be  found  in  [12].  However,  physics-based 
derivation  of  these  coefficients  is  ongoing  task  closely  related  with  numerical  analysis  and  will 
be  presented  in  the  final  report.  The  last  equation  for  has  been  derived  in  [11]  using 

consideration  of  Coleman  and  Noll  on  positiveness  of  entropy  production 

Complete  system  of  governing  equations  needs  to  be  added  to  the  constitutive  relations. 
Combination  of  “elastic”  and  “viscous”  momentum  (or  in  other  words  Newton’s)  equations  for 
each  species  has  the  following  general  form: 


0cr*  d<Ju 


sdiss 


dx^ 


-  +  - 


1 _ =  Y 


DAlamber 


Sx,. 


+  F. 


D'Arcy 


jtD'  Arcy 


V  (-f,  ( A  +  pf'' ))  =  K 

v(-^g  [Pg  +  pT  ))  = 


(27) 


The  second  of  these  equations  we  have  already  used  deriving  novel  generalization  of  D’Arcy  law 
for  unsaturated  porous  media  under  impact  of  dynamic  loading.  Thus,  we  have  described  stress- 
strain  relationship,  forces  causing  these  stresses,  and  constituents  displacements  and  velocities 
resulting  from  these  stresses.  To  complete  this  system  description  one  needs  to  add  the  mass 
transfer  equations  (filtration). 

|-M)  +  V(p,(V,+w,))  =  0 

Finally  we  have  the  system  of  governing  equations  describing  multiphase  flow  in  porous  media 
with  acoustic  streaming.  The  transmission  of  sound  through  unsaturated  porous  media  can  be 
predicted  only  in  the  context  of  the  extended  Biot  model  we  have  developed  here.  The  water  and 
air  move  simultaneously  with  solid  matrix.  The  completed  fully  coupled  model  predicting  waves 
propagations  in  multiphase  porous  media  and  acoustic  streaming  is  given  below: 
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M,-v(v.u,)+m^.v(v.uJ- 


The  next  two  vector  equations  describe  induced  waves  and  acoustic  streaming  in  liquid  and  gas 
phases. 


Close  with  (filtration)  mass  transfer  equations: 


Note  that  when  liquid  and  gas  can  be  considered  incompressible 


As  was  already  mentioned  in  the  section  7  in  accord  with  results  in  hysteresis  in  [12]. 
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There  are  solid-matrix  attenuation  terms,  which  have  exactly  the  same  general  form  as  elastic 
part  but  displacements  are  substituted  with  appropriate  velocities.  As  in  case  of  elastic  solid,  the 
wave  equations  describe  dilatational  and  rotational  waves.  For  the  case  of  strong  coupling 
between  phases  there  exist  two  compression  waves  -  fast  and  slow.  It  is  important  to  note  that 
characteristic  time  scales  for  wave  propagation  and  for  D’ Arcy  mass  transfer  are  very  different, 
which  makes  direct  numerical  analysis  difficult  if  impossible  but  at  the  same  time  allows 
asymptotic  analytical  solution.  What  is  more,  asymptotical  transformations  can  prepare  the 
system  for  numerical  analysis. 

Flow  under  static  pressures  is  expressed  through  “filtration  -  viscous”  velocities  as  follows: 


V  {-s,P,  )  =  {s,y  V,  +  V. 


K. 


gi  ’g 


vK^g) = 


K. 


With  the  constitutive  relations  for  static  pressures  (here,  as  everywhere  in  this  report,  saturation 
S=Ei/<t»): 

P.-P.-J(SVT-^S 

The  general  feature  of  the  developed  model  that  we  account  for  material  morphology  and 
performance  regime  (start  up,  shut  down,  long  steady  working  regime,  etc.)  in  order  to  predict 
fuel  cell  performance  and  degradation  caused  by  increased  water  residual  content. 


9.  Model  summary. 

9.1  Physical  phenomena  included  in  the  model. 

The  model  separates  elastic  and  plastic  properties  of  three  immiscible  phases  involved  in  mass 
and  momentum  transfer  on  the  basis  of  difference  in  time  scale  for  vibration  (elastic  part)  and 
filtration  (plastic  part).  The  plasticity  of  the  solid  porous  phase  is  neglected  as  well  as  all 
temperature  effects. 

The  vibration  part  is  represented  as  three  interconnected  elastic  bodies,  only  solid  having  non¬ 
zero  share  module.  Two  types  of  dissipative  mechanisms  have  been  included:  damping  of  the 
vibration  proportional  to  velocity  gradients  and  momentum  transfer  due  to  phase-to-phase 
interaction. 

The  filtration  part  is  modeled  as  the  motion  in  slowly  varying  static  pressure  fields,  dissipation 
mechanism  being  of  D’ Arcy  type  i.e.  proportional  to  relative  velocities  of  liquid  and  gas. 
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The  above-mentioned  two  types  of  motion  interact  through  dependency  of  saturation  field  at 
every  point.  It  means  that  elastic  parameters  and  slowly  varying  static  pressures  are  functions  of 
saturation.  Two  mass  transfer  equations  for  liquid  and  gas  phases  contain  contributions  from 
both  vibration  and  filtration  parts. 

A  special  equation  that  closes  the  system  Is  a  constitutive  relation  for  two  static  pressures  and  the 
rate  of  saturation  change.  This  equation  determines  the  stationary  state  and,  therefore,  remaining 
liquid  after  all  transient  process  dies  out. 

Of  special  interest  are  non-classical  boundary  conditions  that  result  in  the  net  desaturation  of  the 
porous  sample.  They  should  account  for  dependency  on  the  flux  direction  and  on  the  rate  of  the 
water  removal  from  the  boundary  for  the  net  effect  to  exist. 


9.2  Physical  phenomena  not  included  in  the  model. 

Ideally,  there  should  be  only  one  system  for  momentum  transfer  that  is  not  divided  into  vibration  and 
filtration  parts.  This  could  be  important  in  all  cases  when  time  scales  between  fluids  flow  due  to  external 
pressures  and  caused  by  vibration  are  close  to  each  other.  Such  model  would  require  using  mixed  elasto- 
plastic  constitutive  relations  of  Maxwell- Voight  type. 

Another  potentially  important  gap  to  fill  in  this  model  is  to  evaluate  dependency  of  parameters  in  Leveret 
function  on  vibrations.  That  would  require  deep  analysis  of  droplet  instabilities  trapped  inside  porous 
structure.  The  results  unavoidably  are  morphology  dependent. 


10  Analytical  approach  and  Scale  Separation 


The  primary  goal  of  this  work  is  to  show  that  a  periodic  perturbation  of  the  solid  matrix  or 
pulsation  of  the  gas  phase  causes  non-periodic  flow  of  the  liquid  component.  Although  the 
cornplete  solution  of  the  full  system  of  equations  requires  numerical  simulations  (results  of 
which  will  be  reported  in  the  final  document),  it  is  possible  to  observe  the  effect  using  the  fact 
that  all  filtration  processes  are  substantially  slower  than  the  vibration/pulsation  period.  Under 
such  circumstances  one  can  split  the  complete  system  into  vibration  (Biot)  part  and  filtration  part 

and  solve  them  iteratively.  In  the  first  step  all  saturations  or  in  other  notation 

that  makes  the  whole  system  nonlinear  are 

considered  constant.  The  periodic  solution  of  the  resulting  linear  system  for  |u^,u,,Ug|  can  be 

found.  At  the  next  step,  these  periodic  perturbations  are  imbedded  into  filtration  equations  as  a 
source  term. 


dt 

£ 

dt 


tirP.)+V-(^Pgfg|-(u,-U,)j  =  0 
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We  assume  that  the  periodic  signals  are  in  certain  sense  small  in  order  to  resolve  the  problem 
analytically.  Having  a  small  parameter  at  hand  one  can  use  asymptotic  expansion  based  on 
averaging  principle  [17]  to  obtain  an  averaged  filtration  equations  that  describe  the  slow  non¬ 
periodic  flow  caused  by  periodic  perturbation. 

One  more  important  motivation  for  the  asymptotic  analysis  of  the  derived  equations  comes  from 
the  fact  that  the  accuracy  of  the  numerical  analysis  of  the  wave  processes  is  spatially  limited  by 
the  size  of  the  finite  element  (or  the  step  in  the  finite-difference  scheme).  Therefore,  even 
theoretically,  it  is  practically  impossible  to  develop  a  direct  computational  approach  for  the 
analysis  of  high  frequency  waves  interactions  with  mass  flow  process.  The  only  feasible 
approach  is  to  homogenize  the  mass  transfer  equations  over  high  frequency  acoustic  vibrations 
treated  as  effective  body  forces. 

For  the  sake  of  simplicity  and  without  loss  of  generality  we  assume  that  solid  matrix  is  rigid 
=  const),  both  liquid  and  gas  are  uncompressible  at  the  applied  pressures.  Also,  D’Arcy 
flows  for  both  components  are  uncorrelated.  Crucial  assumption  for  analytical  evaluation,  is  that 
displacements  {u,,U/,uJ  are  caused  by  an  actuator  attached  to  the  porous  matrix  in  such  a  way 

that,  when  it  is  turned  off,  {u,,u,,Ug}  =  0.  So,  the  case  of  small  vibrations,  magnitude  of  the 
displacements  is  proportional  to  small  parameter  |u^,U/,Ug}~A.  Non-destructive  acoustic 
excitation  is  applied  with  A  is  a  small  amplitude  of  external  source. 

Using  equations  (21)  and  derived  generalization  of  the  D’Arcy  law  we  can  write  filtration 
equations  with  periodic  terms  down  as  follows; 
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V  V  J^‘  J 

and  constitutive  relations  for  | : 

£,+e,+e^=\ 

The  following  simplifications  make  the  system  more  tractable  analytically: 

~  lPi>Pt]=  const  -  both  liquid  and  gas  are  uncompressible; 

-  K  J  -  Darcy  flows  for  both  components  are  uncorrelated; 

-  =  const  -  solid  matrix  is  rigid; 

Re-writing  these  relations  through  unique  variable  -  water  saturation  S,  we  immediately  have: 


dt\-S]  e 


where 


//^v.k^(i-5)-’|^v(i-s)[p^-;i.f;]-a.|(i-5).f;‘ 


e  =  \-e^  -  porosity; 


^Is  -V-u^  +  /M,  -V-U/  +£‘m^,  S  V  vt^ 
r.r  •(l-5')-V-U/  +7Wg  -V-Ug 


Here  equation  P  -/»  =  l  +  /rj  has  been  resolved  versus  —  due  to  the 

'  '  5/  0/  j  dt 

simplifications  above,  with  two  branches  correspondent  to  drainage  and  imbibitions  respectively. 
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p.-p,= 


J‘‘’’iS)  +  T-—  for  drainage 
^  '  dt 

f"'”  is) +  T - for  imbibition 

^  ’  dt 


^<0 
\dt  j 

f->ol 

ydt  j 

iK 


As  it  can  be  seen  from  (30),  two  forces  and 


stemmed  from  a  source  of  vibration  affect 


fluid  flow  in  porous  material.  As  they  are  linear  combinations  of  |u^,u,,Ug|  and  their  time  and 

spatial  derivative,  one  can  conclude  that  in  general  the  spectral  composition  of  the  forces 
includes  both  natural  frequencies  of  porous  medium  (Biot  waves)  and  frequencies  that 
correspond  to  the  particular  porous  layer  structure.  The  whole  picture,  therefore,  is  both 
geometry  and  material  dependent  and  requires  numerical  analysis  for  any  specific  geometry.  One 
can  show,  however,  in  very  general  terms,  that  periodic  forces  and  cause  non-periodic 


mass  transport.  Similar  results  have  been  obtained  in  [18]  in  different  setting.  An  approached 
adopted  here  is  that  one  originally  proposed  in  [17]  and  more  generally  in  the  theory  of  normal 
forms  (see  for  example,  [19]).  We  are  looking  for  a  substitution  of  original  variables  such  that  in 
the  new  (slow)  variables  the  mass  flow  equation  does  not  contain  vibration  terms,  both  the 
substitution  and  the  new  mass  equation  being  sought  as  an  asymptotic  expansion  in  the  small 
parameter  A, . 


P,  (x,t)  =  P,{x,t)  +  X-p[  {s,P„J^,x,t)  +  X^-p[  {s,P„J^,x,t)  + ...  (31) 

(x,t) = :^(x,t) + A  • /if  (5,:^,:^,x,t) + • /.f  (^,:^,:^,x,t) + ... 


All  fast  (vibration  terms)  parameters  affect  the  solution  through  the  cross-correlations  averaged 
over  the  time  period.  Substituting  expansions  (31)  into  flow  equations  (30)  and  equating  terms  at 
the  same  powers  of  small  parameter  we  obtain  the  sequence  of  recurrent  differential  equations, 
from  which  all  unknown  functions  of  the  substitution  could  be  defined.  In  other  words,  we  find  a 
solution  for  s^  as  a  functional  of  slow  variables  {S,P,,Pg}.  This  term,  being  purely  vibrational, 

gives  zero  contribution  into  equation  for  the  slow  variables  {S,P,,Pg}.  The  first  non-zero 
contribution  comes  from  ,  that  could  be  found  as  a  functional  of  { S,  i] ,  }  in  the  same  way 

as  s,  .All  details  of  calculations  are  shown  in  the  Appendix  1.  Finally,  the  slow  non-periodic 
flow  is  expressed  through  all  of  these  components  that  play  role  of  body  forces.  Thejnost 
important  results  are  as  follows:  (/)  the  principal  term  of  the  asymptotic  solution  {S,P,,Pg) 

represents  slow  multiphase  flow  under  static  external  pressure  (D’Arcy  flow)  without  acoustic 
impact.  (//)  First  correction  (at  A' )  is  identically  equal  to  zero,  so  there  is  no  first  order  effect  of 
periodic  vibration  on  mass  flow.  In  the  second  -  quadratic  term  the  effect  is  noticeably 
pronounced.  All  vibration  terms  are  coupled  and  their  self  and  cross-correlations  cause 
additional  driving  forces  similar  (in  some  sense)  to  the  effect  of  gravitation  on  the  reservoir  flow. 
The  final  expression  for  slow  non-periodic  flow  is  the  following: 
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dt  Efii 


One  can  see  from  the  expression  above  that  total  non-periodic  flux  is  a  sum  of  the  static 
unsaturated  flow,  vibration  induced  flux  enhancement  due  to  pressure  gradient,  saturation 
gradient  and  their  cross  correlations.  It  is  worth  to  note  that  without  external  pressure  gradient 
the  vibration  does  not  cause  any  non-periodic  effect.  The  third  term  has  always-positive 
coefficient  (permeability  assumed  to  be  a  power  function  of  saturation  [15])  and  increases  with 
inhomogeneity  of  saturation,  i.e.  the  bigger  droplets  the  stronger  effect.  All  terms  include  local 
pressure  gradients,  caused  by  deformation  gradients,  so  the  local  waves-  induced  dilatation 
gradients  are  crucial. 


To  obtain  simple  formulas  for  vibration-induced  contribution  to  the  mass  flow  we  use  method  of 
harmonic  balance.  In  the  simplest  case  of  harmonic  acoustic  excitation,  we  can  express  all 
parameters  as  mono-harmonic  wave  functions’: 

u,  =*  u,  •cos(ft)  r-i-k  x)+’  u,  •sin(<y  r  +  k  x);  i  =  s,l,g; 

and  calculate  additional  pressure  and  permeability  that  this  simple  vibration  induces  in  porous 
material.  First,  a  response  to  this  driving  force  needs  to  be  calculated  in  the  following  form: 

'^.1  ['^.]  [^5,1 

<  p[  >  =  <  ‘'pI  '•cos(<y  /-i-k  x)  +  -  'p[  ►•sin(<y  /  +  k  x) 


Substituting  this  into  differential  ecjuations  obtained  for  and  relating  variables 
{5i(x,r),p'(x,r),pf  (x,/)|,  one  finds  the  linear  system  for  amplitudes 

{  -Sp  Pp  P|,  pf ,  P|*  I  (see  Appendix  2),  from  which  explicit  expressions  for  self  and  cross¬ 
correlations  are  to  be  found.  As  have  been  shown  by  the  method  of  harmonic  balance  the  self¬ 
correlation  of  the  first  asymptotic  term  is  J^^  =||f'  ~  ^u^  All  derived  formulae  are  given 

at  the  end  of  the  Appendix  2  and  the  work  on  the  analysis  is  still  in  progress.  However,  it  is  clear 
that  the  acoustic  streaming  driving  force  is  proportional  to  the  ratio  of  the  fluxes  induced  by 

^  It  can  be  easily  generalized  using  Fourier  series 
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dilatation  and  by  the  saturation  gradient  F  ~  paP' - — ^  +  «k  •  u  .  General  relation  between 

k  •  K;k 

acoustic  streaming  driving  force  and  frequency  depends  on  dependency  between  wave  vector 
and  the  frequency  and,  therefore,  varies  with  geometry.  For  the  general  case  of  the  plates 

dispersion  relation  is  k~  03^  and,  subsequently,  =  f’  -F^  ~  +  ccP^ ,  where 

parameter  5^  is  controlled  by  geometry.  Below  we  show  a  limiting  case  where  an  analytical 
analysis  allows  for  qualitatively  important  results.  Shaking  of  the  structure,  as  a  rigid  body  is  a 
particular  case  we  investigated  during  our  permeability  test  program.  For  this  case  driving  force 
is  proportional  to  the  square  of  the  frequency,  to  the  magnitude  of  the  vibration,  and  to  the  some 
power  of  the  saturation  level:  F  ~  paP  ■  Amp  ■  ,  allowing  estimation  of  the  frequency  role  on 

induced  water  transport. 

We  close  this  section  with  some  concluding  remarks.  Based  on  the  first  results  of  the  asymptotic 
analysis  one  may  see  that  (z)  Acoustic  streaming  effect  in  the  porous  unsaturated  media  does 
exist;  (ii)  Result  depends  both  on  global  geometry  as  well  as  on  local  effects.  The  power  of  the 
source  depends  on  geometrical  parameters,  (in)  Important  [and  now  possible]  to  study  wave 
excitation  for  different  options:  air  pulsation,  structural  waves,  surface  waves,  etc. 


1 1  Experimental  Study 


We  have  conducted  tests  intended  to  measure  liquid  water  permeability  of  the  cathode  layers 
under  vibration,  static  pressure  and  their  combination.  Under  operating  conditions,  water  is 
wicked  away  by  the  pressure  difference  between  catalytic  layer  (where  water  is  generated)  and 
water  transport  plate  (WTP)  with  negative  pressure  (sink).  In  the  development  of  the 

experimental  methodology  we  have  used  and 
I  A  cidir  modified  original  approach  of  UTC  Fuel  Cell  [20]. 


Bolted  down  -  40-60 
psi  equivalent 


Water  in  (Supply 
Figure  8.  Specimen  Setup 


A  scheme  of  the  setup  is  shown  in  Figure  8  and  the 
picture  of  the  stand  is  given  in  Figure  9.  A  sample 
(porous  layer  from  fuel  cell  cathode)  of  two-inch 
diameter  is  placed  between  two  bi-polar  plates.  The 
two  bi-polar  plates  (WTP)  are  then  housed  between 
two  Plexiglas  plates.  0-rings  are  inserted  between  the 
Plexiglas  and  plates  contact  to  provide  a  good  seal. 
One  of  the  plates  is  connected  to  a  water  reservoir. 
The  second  plate  is  channeled  and  the  vacuum  pump 
provides  a  vacuum  (negative  pressure).  The  water 
flow  rate  is  measured  by  the  high  precesion  balance 
with  the  accuracy  of  0.1  microgramm.  We  apply 
vibration  using  shaker  on  top  of  external  constant 
pressure  gradient  to  the  sample.  We  compare  the  flux 
(permeability)  under  different  conditions. 


The  principal  difficulty  of  this  test  is  the  necessity  of  comparison  between  very  small  fluxes  of 
the  order  of  microgramms/sec.  Another  principal  limitation  of  the  developed  setup  originated 
from  small  size  of  the  specimens  (thickness  of  the  porous  carbon  paper  is  about  200  p).  Using 
such  a  test  setup  is  impossible  to  generate  a  dilatational  wave  within  this  porous  structure. 
Taking  all  these  circumstances  into  account  we  have  applied  vibration  (shaking)  of  the  whole 
setup  as  a  rigid  structure.  It  means  that  during  this  type  of  the  tests  we  were  able  to  investigate 
only  the  role  of  inertial  forces  but  dilatation.  Analysis  of  the  role  of  the  dilatational  waves  is  now 
ongoing  task  and  bigger  setup  and  results  will  be  discussed  later.  However,  current  results 
demonstrate  the  non-zero  effect  of  the  zero-mean  vibrations  on  the  permeability. 

We  have  tested  two  types  of  the  porous  substrate  materials  with  different  degrees  of 
hydrophobicity.  The  first  set  of  tests  have  been  conducted  on  10-wt%  PTFE  (10%  of  Teflon),  90- 
wt%  carbon  sublayer.  This  material  manifests  low  hydrophobicity.  We  apply  two  different 
pressure  gradients.  The  first  is  water  pressure  driving  water  into  the  porous  membrane;  the 


Pressure  meters 


Figure  9  Permeability  measurement  test  setup 

second  one  is  the  so-called  vacuum  pressure,  driving  water  out  of  the  porous  plate.  The  porous 
plate  is  not  completely  sealed  in  the  sense  that  gas  pressure  is  assumed  to  be  equal  to 
atmospheric.  Typical  test  results  are  shown  in  Figure  10  when  water  pressure  was  20  Kpa  and 
vacuum  pressure  was  10  Kpa.  During  both  reported  tests,  the  flow  rate  was  steady  and 
approximately  equal  to  0.205  ±0.005  milligrams  per  minute.  The  shaker  were  started  at  the 

point  of  time=4  (red  line  on  the  plot).  Power  of  the  signal  was  500  pV/ms  with  broadband 
frequency  excitation  level. 
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As  one  can  see  from  the  Figure  10  the  increase  of  the  flux  has  been  observed.  If  we  integrate  the 
area  under  the  experimental  curves,  we  obtain  the  total  flux.  Thus,  the  total  flux  increase  due  to 
specimens  shaking  equals  to  0.1  milligrams.  Total  amount  of  additional  released  water  is  about 
15%  of  usual  connected  water  flow.  Estimations  of  the  amount  of  trapped  (disconnected)  water 
based  on  Leverett  function  show  that  15%  is  close  to  the  upper  bound  of  the  possible  trapped 


Figure  10.  Permeability  Measurement  with  water 
pressure=20  kPa;  vacuum  pressure^lO  kPa 

water  amount.  Thus,  during  vibration,  almost  all  trapped  water  has  been  removed.  In  both  tests 
the  the  amount  of  the  released  water  was  the  same.  The  increase  in  water  flux  starts  after  2-3 
minutes  after  turning  the  shaker  on.  Electronic  balance  reading  takes  place  every  minute,  so  the 
difference  is  related  to  the  time-reading  rounding. 

The  second  set  of  test  have  been  performed  on  highly  hydrophobic  materials,  where  50%  of 
teflon  is  added  to  the  carbon  matrix.  The  flux  has  been  unstable  under  water  pressure  of  3  psi 
and  vacuum  pressure  of  1.5  psi.  It  requires  some  time  to  accumulate  water  on  the  specimen 
surface  prior  it  penetrates  into  specimen.  Because  of  high  hydrophobicity,  water  practically  does 
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Figure  11.  Comparison  of  water  flow  in  highly  hydrophobic  porous  matrix  without  (left)  and  with  (right)  vibration 
source. 


not  stop  inside  the  specimen.  As  a  result  one  can  observe  discrete  droplets  coming  out  the 
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hydrophobic  specimen  on  the  graph  in  Figure  11(a).  Initially  the  time  interval  between  droplets 
was  2000  minutes,  and  slowly  decreases  with  each  successive  dropllet  cycle.  It  might  be 
explained  by  the  increased  of  residual  saturation  in  the  porous  media  and  subsequent  decrease  of 
the  Leverett  function.  In  other  words,  with  the  increase  of  the  residual  saturation,  the 
hydrophobic  resistance  to  new  water  imbibition  getting  smaller  and  the  droplets  frequency 
increases.  When  we  turning  the  excitation  source  on,  the  time  interval  practically  immediately 
decreases  and  stabilizes  at  400  minute.  It  demonstrates  the  strong  vibration  effect  on  the  residual 
water  distribution  and  on  the  amount  of  residual  water  saturation  in  highly  hydrophobic 
structures. 

As  a  conclusion  to  this  section,  we  note  that  extremely  precise  experimental  procedure  has  been 
developed  and  conducted.  Test  results  demonstrate  a)  existence  of  the  effect  of  zero-mean 
vibration  on  directional  unsatutrated  water  flow  in  the  porous  structure  (fuel  cell  cathode 
compartment),  b)  instability  in  the  water  flow  in  highly  hydrophobic  structures  under  vibration 
and  the  need  in  wave  generation. 
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12.  Numerical  Implementation  and  Results. 


12.1  Computational  Tool 

At  this  part  of  the  final  report  we  focus  mostly  on  the  Task  3  developing  a  comprehensive 
numerical  model.  We  implement  the  derived  constitutive  equations  into  the  finite  element 
software  FEMLAB  and  conduct  numerical  analysis  with  different  boundary  and  initial 
conditions.  FEMLAB  has  been  chosen  because  of  its  ability  to  combine  physical/mathematical 
problems  arisen  from  different  fields  in  one  generalized  system.  Combination  of  the  different 
types  equations  into  one  template  is  the  major  advance  of  the  FEMLAB.  In  our  problem,  we 
combine  hyperbolic  system  (waves  equations)  for  dynamics  response  with  mass  transfer 
equations.  In  our  system  equations  (26)-(28)  can  be  modified  to  the  system  of  two  elliptic 
equations  and  the  first  order  (transport)  equation  for  water  saturation  S.  We  use  the  advantage  of 
multiphysical  coupling  of  (1)  acoustics  and  (2)  porous  media  multiphase  flow.  FEMLAB  allows 
running  uncoupled  problems  (pure  dynamics)  and  (static  mass-transfer)  automatically  without 
any  reformulation. 

The  underlying  structure  of  FEMLAB  is  a  system  of  PDE,  which  might  be  represented  in  one  of 
the  following  forms:  coefficient  form  (better  suitable  for  almost  linear  problems),  general  form 
( divT  =  f  )  where  matrix  T  and  forces-vector  f  are  explicitly  specified  by  user;  and  weak  form 
that  is  very  convenient  for  coupled  problems.  It  is  important  to  note  that  in  finite  element 
method,  all  equations  eventually  should  be  expressed  in  weak  form  anyway.  Thus,  we  implement 
all  our  equations  as  a  user  formulated  routines  in  FEMLAB  template. 

FEMLAB  integrates  seamlessly  with  MATLAB  that  provides  the  computational  platform. 
FEMLAB  models  run  in  MATLAB  structures  environment  and  we  can  use  all  MATLAB 
programs  and  procedures.  As  our  numerical  experimentation  has  demonstrated,  the  major 
disadvantage  of  the  FEMLAB  is  underdeveloped  meshing  procedure  and  impossibility  to  modify 
the  solver  routine.  It  limits  the  applicability  of  this  tool  (in  it  current  version)  to  relatively  simple 
geometries  with  not  large  number  of  elements. 


12.2  Formulation  and  Normalization  of  the  Governing  Equations. 

First,  we  consider  the  filtration  part  of  the  equations  as  it  follows  from  Eq.  (21) 

^  +  V-(w,)  =  0;  and  w,  =5-(v, -vj 
ot 

This  is  expression  for  water  content,  similar  equation  can  be  easily  written  for  gas  phase. 

For  the  successful  numerical  analysis  we  have  to  determine  material  constants  and  extract  J- 
Leverett  functions  from  the  available  experimental  and  theoretical  evaluation.  We  conduct 
combined  theoretical-experimental  study  on  formulating  and  calibrating  models  parameters.  The 
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particular  form  of  J-Leverett  function  and  its  relation  with  acoustic  loading  is  the  subject  of  the 
thermodynamic  analysis. 

Now,  have  to  emphasize  the  very  important  point  of  the  “elastic-plastic”  decomposition.  This 
has  been  already  described  in  the  “theoretical”  part  of  the  report  but  we  repeat  it  again  in 
numerical  part  of  the  report.  Motion  of  the  fluid  is  to  be  decomposed  into  two  parts:  elastic  or 
dynamically  induced  inertial  motion  and  the  motion  induced  by  pressure/saturation  and 
controlled  by  modified  D’Arcy  law.  In  other  words,  we  decompose  the  total  motion  into 
“elastic”  wave-related  motion  and  “plastic”  filtration  flow*: 

^  ^d}'namicat  filtration 

The  D’Arcy  part  of  the  flux  is  described  by  gradients  and  as  shown  in  Chapter  7  can  be 
expressed  as 


+  V-(5-(v,-vJ) 


^D'Arcy-  ^  ~  >  ^ - 

S  /// 

if  we  neglect  effects  of  solid  finite  dilatation.  Finally,  the  conservation  law  can  be  re-written  as 
follows: 

f  =  v{*(5)^].V,5.(v,-v.), 

The  similar  equation  can  be  derived  for  the  gas  phase.  Equation  for  the  dynamic  equilibrium  .has 
the  following  general  form 

r(5)~-/>«+/>'=/(5) 

ot 

These  constitutive  relations  are  supplemented  by  the  definition  of  the  Leverett  function, 
simplified  form  of  which  is:  J(5)  =  — ^ ,  where  S,eft  is  the  left  saturation  asymptote 

and  Sin  is  the  dynamic  saturation  equilibrium  specific  for  the  pair  porous  media-fluid.  In  our 
current  calculations  T(S)  has  been  taken  as  a  constant.  In  limiting  case  of  equilibrium  (—  =  0) 
this  equation  is  equivalent  to  well-known  Laplace  —Leverett  equilibrium  condition. 

Combining  together  both  transport  equations  for  liquid  and  gas  and  the  equilibrium  condition 
and  eliminating  time  derivatives  from  the  first  two  equations,  we  obtain  the  following  system  of 
two  elliptic  and  one  first-order  equation 

-T{S)K,V-[s^.VP,+S-VSP,]  +  VS(v,-vJ  +  P,-P,=-J(S) 
--r(5)-/:,.V.[(l-5)^V/>-(l-5).V5.P,]  +  V.(l-5)(v^-v,)  +  />,-i>=J(5) 
Pt-P2+J(S)  =  T(S)—;  =  — 


This  decomposition  is  also  correct  for  finite  deformation  and  similar  to  the  velocity  gradient  decomposition  in 
deformation  plasticity. 
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This  system  describes  only  “viscous”  part  of  the  system,  or,  in  other  words,  only  filtration  part.  It 
is  coupled  with  dynamic  part  mostly  through  species  velocities  but  also  through  the  pressure 
distribution. 


12.3  Normalization  and  Calibration  of  the  Parameters 

During  this  project  we  analyze  harmonic  signals,  so  the  natural  choice  for  time  scaling  is  r  = 
For  the  stability  of  numerical  procedure,  it  is  desirable  to  get  the  coefficient  in  front  of  highest 

X  X 

derivative  of  the  order  of  unity.  This  dictates  the  choice  of  special  scaling:  x  =  -j=  =  . — — . 


This  set  of  substitutions  leads  to  introducing  of  new  normalized  parameters/variables: 

C  P 

C  =  and  P  =  — .  It  is  easy  to  see  from  these  groups  that  the  role  of  applied  “static” 

<oT  oT  .  .  9 

pressure  gradient  decreases  with  the  increase  of  the  applied  acoustic  frequency  . 

There  are  two  limited  cases  of  the  problem  of  acoustic  enhancement  of  mass  transfer  in  porous 
media.  The  first  reflects  the  case  when  small  vibration  perturbs  filtration  process.  The  principal 

term  in  the  equation  here  is  V  •  [  k{S)  .  This  corresponds  to  the  numerical  situation 

when  coefficient  of  filtration  K  is  relatively  large.  Asymptotical  analysis  of  this  situation  has 
been  reported  in  the  section  above.  The  other  limiting  case  corresponds  to  small  permeability 
(say,  A:~le-15)  and  the  problem  becomes  close  to  singular  perturbation  (fx  +  x  +  a>^x  =  / ).  The 
numerical  finite  element  software  developed  in  this  program  covers  wide  range  of  the  parameters 
variation. 

In  order  to  calibrate  and  verify  the  developed  model  we  need  evaluate  material  parameters  and 
fit  model  parameters  against  experimental  data. 

For  isotropic  Toray  paper  the  following  material  parameters  are  accepted: 


X  =  «  30  Ksi\  n  =  «  IQKsi-,  p  « 1000%/ and  porosity  (/> »  0.75 . 

Material  properties  for  Toray-paper  bulk  compressibility  is  »  30  Ksi  and  for  structural  fibers 
is  »  20  Msi ,  also  »  25  Msi « 1 60  GPa .  Viscous  properties  and  acoustic  waves  decay 

for  this  material  has  been  accepted  from  literature  review  and  based  on  averaging  for  test  results 
for  different  woven  carbon  materials: 


^  Other  normalizations  are  possible,  however  we  stopped  on  the  described  above.  We  have  checked  several  different 
substitutions  of  which  the  most  numerically-convenient  lead  to  the  normalization  of  the  coefficient  in  front  of 

(X  +  2p)  k,  __x  , 

D’Arcy  term  in  the  dynamic  system.  For  example,  r  = - - j - jyi;  where  /  is  the 

/  Pi  / 

characteristic  size 
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wave  decay  (y)  based  on  Delany  and  Bazley  laws  we  were  able  to  estimate 

T] «  0.05  ^  0.5  Mpa  •  s  and  ^  «  0.02  0.2  MPa  •  s .  Elastic  generalized  Biot  coefficients  have  been 

evaluated  based  on  elastic  properties  as  follows: 

q  =  690  MPa  •  t  =  0.04  MPa  -(f)-,  r  =  220  MPa  •  5  =  0.12  MPa  ■  (f>\  where  modified  Biot’s 

coefficients  has  been  expressed: 

Q-q£i  =  q(f)S',  R  =  r£,  =  r^S;  T  =  t£^  S);  5®"”  =  s£^  =  5^(1  -  S).  It  is  important 

that  we  assumed  that  q,  t,  r,  and  5  are  constants  and  do  not  depend  on  saturation. 

Next,  we  estimated  transport  properties,  which  requires  fitting  of  model  parameters. 

Viscosity  of  water  is  assumed  to  be  equal 


=  O.Scp  =  5e  -  4  ^  ^ =  0.022cp  =  2e  -  5  ^ .  We  need  to  evaluate 

permeability  of  water  in  the  porous.  Using  Klinkenberg  effect  we  approximate  ~  1.4  • 

Testing  of  water  permeation  through  substrate  at  different  pressure  gradients  have  been 
conducted  at  UTRC  and  UTC  Fuel  Cell.  We  fitted  the  flux  vs.  pressure  finite  element 
calculations  against  experimental  data  as  shown  in  Figure  12  below. 


Figure  12.  Fitting  transport  parameters  against  experimental 


We  varied  several  major  parameters  of  the  system  aiming  to  obtain  the  closest  fit.  Leaving  the 
question  of  the  uniqueness  open,  we  are  able  to  find  the  set  of  parameters  providing  good  fit. 
With  decrease  of  viscosity  T,  the  flux  increases;  Coefficient  Co  -magnitude  of  Leverett  function 
defines  the  rate  of  stabilization  and  with  the  increase  of  Co  the  flux  decreases. 

Calculated  fitting  parameters  are  shown  in  Table  1. 
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More  exact  fitting,  especially  for  big  pressure  gradients  can  be  obtained  by  the  use  of 
Forchheimer  model'®.  We  have  done  some  numerical  experiments  in  that  direction  but  this  work 
is  out  of  scope  of  the  reported  project. 


12.4  Boundary  conditions. 

The  system  derived  above  belongs  to  mixed  hyperbolic-parabolic  type,  vibration  part  being 
hyperbolic  and  filtration  parabolic  wit  ordinary  differential  equation  for  static  pressures. 

For  vibration  part  zero  normal  forces  are  defined  on  free  boundaries: 
rij -a;  =(X-0,+M,, -e,  -e,.  =0 

"y  =~Pg  ‘"y  ="y  '^v  +^gi  ■^1)  =  ^ 

on  fixed  boundaries  zero  displacements  are  defined  for  all  phases: 
u,  =0;u,  =0;u^  =0. 

Boundary  with  excitation: 

(n-Uj)n  =  f(/);u^  -(n  •u^)n  =  0;u,  =0;u^  =0 

where  f(/)  is  a  given  motion  of  the  excitation  boundary.  It  can  be  pulse-like  or  periodic. 

To  be  consistent  it  requires  appropriate  definition  of  velocities: 

(n-vjn  =  5f(0/5t;v, -(n-vjn  =  0;v,  =0;v^  =0. 

For  filtration  part  the  boundary  conditions  should  reflect  the  physical  conditions  of  water 
removal  from  the  boundary.  As  one  can  see  from  the  Constitutive  relations  for  static  pressures 
above,  the  sign  of  derivative  dSjdt  on  the  boundary  depends  on  the  relation  between  boundary 

values  of  static  pressures  ,  which  means  that  varying  due  to  vibration  saturation  S  might 

change  this  condition  from  imbibitions  to  drainage  and  vice  versa.  To  obtain  net  drying  effect  a 
boundary  condition  that  prevents  liquid  flow  from  the  outside  must  be  imposed  when  dSj dt  >0. 
It  also  means  that  to  achieve  overall  drying  effect  caused  by  vibration,  a  physical  mechanism 
preventing  reverse  flow  back  into  sample  must  be  found.  In  other  words,  a  boundary  that  is  open 
for  liquid  exchange  must  favor  liquid  flow  in  one  direction.  That  can  be  due  to  fast  removal  of 


The  transport  law  for  such  a  model  has  the  following  form: 


= - predicts  the  convex  shape  of  the  flux  vs. 

V  /  I  j 


pressure  curve. 
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water  from  the  boundary  with  the  rate  higher  than  the  frequency  of  vibration  or  applying  special 
hydrophobic  boundaiy  layer.  Therefore,  the  boundary  conditions  on  the  boundaries  opened  for 
liquid  flow  is: 


P  =  P  ■  P  =  p 

g  boundary  ^g^^l  boundary  *l 


//|5<0 


J,;if  0 


boimdary  dg  -  J,  -  0  moaus  the  best  case  condition  for  acoustic  water  removal. 


^  Flux, 

11 

where 

1  1 

1 2.5  Results  of  Numerical  Calculations  for  Different  Boundary  Conditions 

First  we  considered  the  problem  in  the  rectangular  with  pressure  defined  as  boundary 
conditions  as  shown  in  the  scheme.  Vibration  is  applied  to  the  upper  boundary  as  follows: 

=  Asm{o)t-kx),  and  .  There  are  two  limiting  cases  depending  on  material 


Figure  13.  Deformed  shape  of  the  plate  . 
Contours  of  Saturation  are  shown. 


parameters  and  BC.  The  first,  when  vibration  [acoustic  pressure]  is  small  with  respect  to  applied 
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to  boundaries  “static”  pressure  gradient.  The  second  case  corresponds  to  the  situation  when 
acoustic  pressure  is  large  in  comparison  with  static  pressure  gradient.  We  have  considered  both 
situations.  The  Figure  below  illustrated  the  simplest  deformed  shape  of  the  structure. 

The  results  of  calculations  corresponding  to  the  first  situation  are  shown  in  Figures  13-16. 
Saturation  distributions  under  both  static  and  small  vibration  conditions  are  similar  to  each  other 
and  show  steep  saturation  gradient  near  to  sink.  Figure  14  illustrates  formation  of  steady 
distribution. 


Figure  14.  Saturation  Distribution  in  the 
layer  after  different  time  till  reaching  the 
steady  regime.  No  vibration  applied. 


Figure  15.  Distribution  of  Steady 
Saturation  with  and  without  vibration 


The  total  effect  of  small  vibrations  on  total  water  content  is  small  and  illustrated  by  graphs  in 
Figure  15.  These  direct  finite  element  calculations  are  in  good  agreement  with  analytical 
(asymptotical)  analysis  (part  9),  which  in  turn  predicts  second  order  effect  of  small  vibration  on 
water  content  in  the  porous  unsaturated  media.  It  is  important  that  in  this  calculation  the  water 
transport  is  primary  controlled  by  filtration  process. 


The  most  pronounced  effect  of  vibrations  on 
transport  process  is  substantial  acceleration  of 
water  transport.  Under  vibration  impact, 
transient  period  is  insignificant  in  comparison  to 
pure  D’Arcy  filtration  (static)  case.  Thus,  under 
pressure  boundary  conditions  and  dominating 
filtration  term  the  water  just  follows  the  motion 
of  the  solid  phase.  This  leads  to  fast  process 
stabilization  and  to  small  integral  effect  in  flux 
evaluation  as  shown  in  Fig.  16.  The 
homogenization  of  water  distribution  is 
significantly  faster  under  vibration  impact. 
Increase  of  the  effect  relates  with  the  increase  of 
the  term  V  •  (5(v„„„,  -  J) .  Therefore,  this 


Figurel6.  Decrease  of  transient  time 
with  frequency  of  applied  vibration 


product  after  averaging  should  be  not  only  significant  but  noticeably  changes  with  space  to 
guarantee  existence  of  not  small  derivative  (div). 
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MacOJii 


Figure  1 7.  Formation  of  the  boundary  layer  on  the 
top  boundary  under  acoustic  impact. 


During  each  half-cycle  water  efficiently  drains  from  the  bottom  boundary  and  imbibes  back 
during  another  half  cycle.  The  same  situation  is  observed  on  upper  boundary.  Only  with  the 
development  of  the  mechanism  preventing  back  imbibing  from  the  surface  or  in  other  words  the 
way  of  fast  removing  water  from  the  surface  would  lead  to  effective  acoustic  drying. 

Another  qualitative  situation  appears  when  filtration  coefficient  is  small  (as  typical  for 
Toray  paper  and  bi-layer)  and,  therefore,  coefficient  in  front  of  highest  derivative  is  a  small 
parameter.  This  and  pressure  boundary  conditions  leads  to  creation  of  a  boundary  layer.  The 
characteristic  size  of  such  boundary  layer  can  be  evaluated  by  equating  orders  of  both  terms  in 

k 

transport  equation  above:  I  ~  KT  ~  ~  The  thickness  of  the  boundary  layer  is 

^water 

about  several  percents  of  the  cathode  porous  layer.  All  changes  consisting  the  boundary 
conditions  with  the  bulk  solution  take  place  in  this  boundary  layer  and  saturation  in  the  internal 
part  of  the  porous  plate  practically  does  not  change  with  vibrations.  As  can  be  seen  from  the 
Figure  17,  the  region  with  decreased  water  content  occurs  inside  the  porous  plate  and.  the 
saturation  is  higher  at  surface  regions,  therefore  water  flux  is  directed  inside  the  porous  plate.  In 
other  words,  during  vibration,  plate  starts  imbibe  water,  if  pressure  is  defined  on  the  boundary. 
This  is  very  important  situation  for  PEM  FC  materials  because  calculations  fit  range  of  material 
constants.  This  leads  to  the  conclusion  that  the  most  important  process  of  acoustic  drying  is 
effective  water  removal  from  the  surface.  Assuming  it  is  technically  possible  to  realize,  we 
analyze  water  transfer  under  this  semi-permeable  boundary  conditions.  The  modeling  results 
provide  us  with  upper  bound  estimations  of  the  amount  of  acoustically  induced  water  flux. 
Acoustically  induced  flux  is  naturally  more  pronounced  if  applied  static  pressure  gradient  is 
small.  We  have  numerically  analyzed  in  details  the  situation  when  the  in-bound  water  flux  at  the 
sink  boundary  is  semi-sealed.  As  was  described  in  the  previous  section,  we  implemented  such  a 
semi-permeable  physical  model  by  constraining  the  boundary  conditions  as  j  TndA  >  0 . 

boundary 


The  set  of  cross-section  snapshots  of  water  saturation  (S)  reflecting  dynamics  of  saturation  with 
normalized  time  under  influence  of  structural  waves.  “Low”  bottom  is  semi-sealed:  if 
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f  ^^>0,  then  flux  is  zero,  otherwise  is  not  constrained.  Pressure  on  the  ^^upper” 
J  dv 

boundary 

boundary  has  been  set  constant  P=0.  The  detailed  (in  time)  set  of  snapshots  is  shown  in  the 
Appendix  III.  Below  we  show  only  the  most  important  pictures. 


Figure  18.  Series  of  snapshots  of  water  saturation  distributions  under  the  impact  of  structural 
vibrations.  Development  of  thin  “boundary  layer  on  the  top  boundaries  observed  under 
“tension”  mode  of  vibration. 


As  can  be  seen  from  the  Figures  18  reflecting  the  acoustically  driven  mass  transfer,  the 
boundary  layer  is  developed  on  the  boundary  where  the  pressure  has  been  defined.  Also,  the 

integrated  amount  of  water  in  the 
plate  gradually  decreases  on 
approximately  5%  and  stabilizes 
there.  Process  starts  with  relatively 
fast  drainage,  then  after  first  half  of 
vibration  period  the  internal  pressure 
gradient  is  developed  in  such  a  way 
that  water  would  be  imbibed  from  the 
upper  boundary.  Up  to  time  ~  4.8 
semi-permeable  boundary  conditions 
lead  to  invariable  amount  of  water 
content  in  the  porous  media  (see  Fig 
19).  Then,  during  another  period  of 
vibrations,  pressure  conditions  change 
again  and  drainage  takes  place.  After 
this  the  process  is  almost  stabilized 
and  further  drainage  has  not  been 

numerically  observed.  The  most  important  conclusion  from  these  calculations  can  be 
summarized  as  follows; 

(i)  The  acoustically  driven  water  removal  from  the  porous  plate  is  numerically  observed; 


Figure  19.  Dynamics  of  acoustically  driven  water 
drainage  through  semi-permeable  boundary 
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(ii)  Acoustic  drying  is  possible  if  and  only  if  water  removal  from  the  surface  mechanism 
does  exist.  Otherwise,  water  is  imbibed  back  immediately. 

(iii)  The  acoustic  drying/drainage  is  very  fast  process,  which  requires  only  several 
vibration  periods  to  reach  saturation; 

(iv)  Total  amount  of  the  removed  water  under  given  parameters  is  not  exceed  5%' ' . 

From  the  series  of  snap-shots  showing  changes  in  liquid  and  gas  pressure  with  structural 
(solid)  vibrations  demonstrate  the  structure  of  acoustic  waves  in  the  liquid  and  gas  phases  of 
partially  saturated  porous  media.  Model  calculations  shown  in  Figures  20  (Detailed  series  of 
snap  shots  are  given  in  Appendix  III)  clearly  show  that  liquid  and  gas  motions  follow  the 
solid  structure  vibrations.  This  effect  takes  place  solely  due  to  the  “effective  viscosity” 
generated  by  D’Arcy  term  in  the  constitutive  equations. 


Figure  20.  Series  of  snap-shots  demonstrating  changes  in  liquid  and 
gas  pressure  with  structural  (solid)  vibrations 


Figure  21  illustrates  in-phase  vibration  of  liquid  and  solid  phases  in  both  structural  vibration 
modes:  longitude  and  shear  waves.  Therefore,  the  dynamics  of  saturation  (liquid  and  gas)  within 


"Parameters  variations  might  lead  to  increase  of  the  removed  water,  but  it  would  not  exceed  10%. 
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the  porous  structure  is  a  strong  function  of  structural  waves.  Very  important  result  that  can  be 
deduced  from  this  numerical  modeling  is  that  under  structural  vibrations/waves  applied  to  the 


Figure  21.  Displacement  along  (a-left  four  figures)  x-direction  (u,)  of  solid  (left)  and  liquid  (right)  phases  and 
(b-  right  four  figures)  y-direction  (U2)  of  solid  (left)  and  liquid  (right)  phases 

partially  saturated  porous  body,  the  liquid  and  gas  pressures  have  a  non-trivial  form  and  shape 
within  the  body.  The  pressure  strongly  affects  the  diffusion  coefficients  and,  therefore,  affects 
transport  of  species.  The  gas  (oxygen  in  our  case)  transport  would  be  determined  by  not  a 
diffusion  as  previously  been  accepted  but  by  convection-diffusion  mechanism. 


13.  Concluding  Remarks 

We  have  developed  a  unified  model  of  structural/acoustic  wave  propagation  in  the  PEM  cathode 
compartment  and  coupled  it  with  mass  transfer  in  the  porous  media.  This  is  the  first  known  to  the 
authors  model  of  such  kind  coupling  acoustics  and  filtration  in  porous  media.  A  set  of  consistent 
constitutive  relations  and  governing  equations  for  multiphase  partially  saturated  porous  dynamic 
media  have  been  developed  using  Lagrangian  approach  with  subsequent  inclusion  of  various 
dissipative  mechanisms.  Developed  unified  methodology  is  useful  for  the  wide  range  of  Fuel 
Cell  problems  as  well  as  for  the  wide  range  of  other  porous  structures  and  could  serve  as  an 
important  design  tool. 

It  has  been  demonstrated  that  the  phase  saturations  have  strong  impact  on  wave  dynamics  in 
porous  continuum.  Explicit  expressions  for  generalized  multiphase  Biot-type  coefficients  that 
explicitly  depend  on  the  phase  saturations  have  been  obtained. 

We  have  succeeded  to  generalize  filtration  (D’Arcy)  law  on  dynamic  system  with  water 
saturation  varying  with  space  and  time.  It  has  been  shown  that  the  viration-induced  flux  is  a 
strong  function  of  solid  matrix  distortion,  saturation  gradient  and  inertial  forces.  We  also  were 
able  to  connect  hysteretic  phenomena  to  relaxation  process  due  to  gradients  of  filtration  rate. 

The  developed  homogenization  technique  is  crucial  for  numerical  analysis  of  high  frequency 
impacts.  It  allows  circumventing  restrictions  on  maximum  allowed  number  of  finite  elements  to 
resolve  spatial  frequency  of  the  acoustic/structural  wave.  Using  asymptotic  analysis  we 
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demonstrated  the  crucial  role  of  cathode  geometry  on  the  water-transport  enhancement  effect. 

We  also  obtained  analytical  results  evaluating  the  role  of  acoustic  loading  on  the  flow  and 
frequency  and  wave  number  dependencies  of  the  effective  mass  driving  forces.  It  has  been 
demonstrated  using  averaging  principle  that  vibration  gives  rise  to  net  change  of  saturation 
inside  porous  medium.  This  effect  is  analogous  to  the  known  phenomenon  of  acoustic  streaming 
[18, 21]  but  has  not  been  discovered  before. 

We  conduct  the  extremely  precise  experimental  procedure  to  measure  bi-layer/substrate 
permeability  under  vibration.  Test  results  for  different  materials  demonstrate  (i)  existence  of  the 
effect,  (ii)  instability  in  highly  hydrophobic  structures  under  vibration  and  the  strong  correlation 
between  induced  mass-transport  and  spatial  inhomogeneity  generated  by  structural  waves. 

The  most  useful  application  of  the  equations  derived  for  time-averaged  fluxes  is  when  the 
leading  zero-order  approximation  is  vanished.  In  the  opposite  case,  vibration-induced  fluxes  are 
small  and  the  transport  is  determined  by  static  pressure  gradients.  This  also  specifies  potential 
application  of  the  vibration-induced  saturation  dynamics.  It  is  most  useful  when  for  some  reason 
pressure-induced  transport  is  either  impossible  or  undesirable.  Numerical  results  demonstrate 
that  applied  structural  vibration/acoustic  loading  (i)  intensifies  the  process  dynamics  or  in  other 
words  decreases  time  of  transient  process;  (//)  might  drain  out  up  to  5-10%  of  distributed  in 
porous  media  water  if  the  mechanism  preventing  reverse  acoustic-driven  imbibitions  is  in  place; 
{Hi)  The  fast  water  removal  from  the  surface  is  crucial  for  acoustic-driven  drying;  (/v)  In  a 
working  Fuel  Cell,  vibrations  cause  increase  of  the  water  flux  from  the  porous  cathode  in  both 
directions  out  to  WTP  and  back  to  the  catalytic  layer,  which  in  turn  makes  the  Fuel  Cell 
performance  worse;  (v)  Acoustic/Structural  vibrations  accelerates  the  stabilization  of  the  system 
(porous  media-liquid-gas)  at  some  level  predefined  by  material  properties  in  general  and  by 
Leverett  function  in  particular.  Therefore,  for  over-saturated  porous  media,  vibrations  lead  to  fast 
drying  up  to  stable  level  and,  subsequently,  for  FC  performance  improvement.  If  the  initial 
saturation  is  smaller  then  the  steady-stable  level,  the  porous  layer  would  be  immediately  imbibed 
up  to  the  stable  level  and,  subsequently,  performance  would  deteriorates.  Both  types  of 
phenomena  have  been  observed  (see  also,  [22])  in  practical  nature. 

On  the  base  of  the  results  mentioned  above,  we  make  the  following  practical  conclusions  focused 
on  results  industrial  (fuel  cells)  implementation:  (i)  Acoustic/structural  vibration  drying  can  be 
used  when  other  mass-transfer  driven  processes  are  small.  It  means  that  the  best  operational 
moment  for  the  application  of  the  acoustics  is  the  shutdown  drying  preventing  from  possible 
frosting  of  the  PEM  fuel  cell  cathode  compartment.  Another  optimum  performance  regime  is 
fuel  cell  dry  startup,  when  vibrations  would  significantly  cut  the  start-up  stabilization  (warming) 
time,  (ii)  Fast  and  efficient  water  removal  from  the  surfaces  and  interfaces  is  crucial  for  the 
acoustic  drying.  Much  work  needs  to  be  done  to  develop  such  a  material/technology  design.  This 
work  shows  the  direction  for  the  optimal  material  selection,  (iii)  Based  on  numerical  results,  one 
can  conclude  that  the  acoustically  driven  material  drying  (water  drainage)  is  the  fast  procesi 
therefore,  the  application  of  vibration  field  should  be  performed  as  a  short  impulse  series, 
preferably  while  source  or  sink  of  water  is  sealed  (start-up  and  shutdown),  (iv)  Project  results 
show  that  cathode  compartment  acoustic  cleaning  would  be  beneficial  if  the  cell  is  tumed-off  at 
the  moment  of  cleaning.  Therefore  subsequent  short-time  (~1  min)  disconnections  of  one  or 
several  cells  from  the  stack  for  a  cleaning  would  provide  positive  performance  effect. 
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Appendix  I 

Introducing  the  new  variable  S'  =  £■,  /(l  -  the  system  takes  the  form: 


dt  i^-S  J  e 

where 


•V.K,S-MVS[/>-A-i<'’]  +  A~S.F,^ 


;.;‘.V.K,(1-S)-*  v(i-s)[i>,-;i-F/]+2.-(i-s)-F; 


s  =  \-E^  -porosity; 


f  g-V-u^ +r-V-u, +s-b-S -V -Ug 
F/  r[t-y-u^+s-b-{l-S)-V-u,+s-V-u^ 


I”,  A  +«  ((i-s)" 

and  constitutive  relations  for  | : 

f)  f  3  ^ 

P  -  F  =  (S)  +  r — S  for  drainage  — S  <  0 

^  ‘  ^  ’  dt  [8t  J 


P  -P.  =J‘'”^  (S)  +  T-—S  for  imbibition  [—S  >  o| 

^  ‘  ^  ’  dt  \dt  ) 

f  Q  Q  \  Q 

Here  equation  P,-Pi=Jisi,£o]  +  F  has  been  resolved  versus  —  S  due  to  the 

®  V  */  dt  J  ot 

simplifications  above,  with  two  branches  correspondent  to  drainage  and  imbibitions  respectively. 

Additional,  substantial  for  analytical  evaluation  assumption,  is  that  the  vibrations  {u^,U/,Ug}  are 

caused  by  a  vibrator  attached  to  the  porous  matrix  such  that,  when  it  is  turned  off, 
{u,,u,,uj  =  0.  So,  the  case  of  small  vibrations  {u,,u,,uj~  iis  considered  where  A  is  small 

amplitude  of  external  source. 


Asymptotic  analysis. 
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As  it  can  be  seen  from  the  above,  two  forces  F,Y  stemmed  from  a  source  of  vibration  affect 
fluid  flow  in  porous  material.  As  they  are  linear  combinations  of  {u^,u,,Uj}  and  their  time  and 

spatial  derivative,  one  can  conclude  that  in  general  the  spectral  composition  of  the  forces 
includes  both  natural  frequencies  of  porous  medium  (Biot  waves)  and  frequencies  that 
correspond  to  the  particular  porous  layer  structure.  The  whole  picture,  therefore,  is  both 
geometry  and  material  dependent  and  requires  numerical  analysis  for  any  specific  geometry.  One 
can  show,  however,  in  very  general  terms  that  periodic  forces  F,Y  cause  non-periodic  mass 
transport. 

A  change  of  variables: 

(5.^, /)+... 

(«.')  =  ^  •  P,’  +  X'  ■  Pi  + ... 

is  sought  in  form  of  asymptotic  expansion  in  A.  such  that 


1  fJ- 


and  |S'(x,/),f’  (x,/),P^  satisfy  filtration  equation  without  vibration  terms: 


S,P„P^,\,t^ 

p' 

's,P„Pg,x,t) 

pf 

5  5  K. 


5/  -5  e 


[k*  (5,^,^)  +  A .  Kf  (5,^,^)  +  Kf  (5,:^,^)...] . 

('-«)"  (v  ■('  -  s)[K  (5.^ .^ )  *X-P‘{U,J.]^^X‘-  P’  (5,^,^)...]) 


Pg  Pi  -  {s,P,,P^^-^-  A-  jf  {s,Pi,P^-k- ■J2{s,P,,P^  +  T~S  fordrainage  f^5'<0 


+  +  +  for  imbibition  (-^5 >0 

'  '  ^  '  dt  {dt 
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The  resulted  equation  for  non-periodic  variables  |5'(x,t),P,(x,t),Pg  (x,t)|  has  the  same  type 
and  shape  as  the  equation  when  vibration  is  absent  =  0) .  It  is  therefore,  possible  to  interpret 
as  an  effective  increase  in  permeability  and  Pl'^{s,P„Pg^  as  an  effective 
increase  in  hydrostatic  pressure  due  to  applied  vibration. 

Functions  K'’«  and  can  be  determined  by  substituting  {5(x,t),i’ (x,t),Pg  (x,t)}with  their 
asymptotic  expressions  through  |s(x,t),i^(x,/),Pg(x,t)J  and  (x,/),;?' (x,t),/?f  (x,/)|  and 
applying  conditions  that  the  time  and  space  averages  of  {5',  (x,t),/7'(x,t),/7f  (x,r)|are  equal 
zero; 


57 


[-S  -  A  •  5,  (5,:^,:^,  X,/)  -  A'  •  5,  (s,;^,  j^,  x,/)  -  ...I  “  ^ 


_ 


[K;(s,i:,7>,)+A-K;(s,7>,fj+2>-K;(s,/;,/>j...]. 
s'  (v  .s[K  (s,^ ,^)+ A .  p;  (S,:^,'^) + A> .  ^  (s,fl,^)...]) 


('[K*(s,/i,/>,)+A.Kf(s,;>,p^)+A’-Kf(s,/;,fJ...].  'I 

^iS  +  A  •  5|  +  A^  •  ^2  •••Jj  • 

///■'•V-  ^  ,fv.(5  +  A.5,+A^52...)[P/+A.p,'  +  A^p^..-A.(f;'’+A.'F/)l 

(s  +  A-5,+A'-J2...)  p  ^  U 


"k,(1-5-A-5,-A^.52...);- 


/^;'v- 


(1-5-A.5,-A'.52...) 


•.2...)  {-HK^^-'K)  J 


_a._(i-5-a..,..)(f-+a.'f;<) 


where 


P,'’'\  (  ^-V-u  +r-V  u, +£'-fc-5-V-u 

F/J  [/.V.u,+^.fe.(l-s).V.u,+5.V.u  I’  [f^ 
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r  q-y-\  +  r-V-\  +  s-b-S-^-'u^ 
K  J"l-f-6-5,-V-uJ'^L-V-'u,+.e-^)-(l-5)-V-y +5-V-’u^ 


I'lpK'l  -/’A  |rf  (“/-“.)  f  P,|'y+(«A  +  P,('S''-i)-^«<)|-{'"/-'".) 

-“.)^  [p.  I'u,  +(p.  A.  ((i-s)"' -i)- A,)|('«, -'», 

(Here  j’u./u^'ujare  the  second  order  terms  in  X  to  be  calculated  from  the  vibration  part  of 
the  foil  equation  set  with  substituted  with  je,  +  £■  •  5,,^^  -  £•  ■  -  not  needed  in  equations 

for  slow  variable  in  the  second  order  in  X ). 

'P,-li+X\p^-p[)  +  X^\pl-p[)  = 

(s)  .  A . + -yi..  j  +  7- ■  +  A .  +  A  .  .4  = 

Jf  (s)  +  A  •  Jf  (5,^,^)  +  X  ■  Jf  (5,:^,^)  +  T  •  ■|:(5  +  X-s,+X^-s^)  for  drainage  j^— 5  <  0^ 


Pg-Pi+x-  [pf  -p[)+x'^-  [pI  -  y )  = 


a/""* (5)  5V"”'’  5J"”*(5) 

W  QS  55  2  55 


5V"”"  5  (e'j^  5J"”"  5  Q,-  ^  . 

2.  - + - X-Ls,  +T-^S  +  X-s,  +  X^-sX 

2  55  ^  dF  ' 


f"’’ (S^  +  X-  Ji"’’’ (s,P„pF)  +  X  -  Jf  (s,Pt,PF)  +  T -^^(S  +  X- s^  +  X'  ■  for  imbibition  >  0 


Now,  collecting  terms  with  the  same  power  of  X" ,  one  gets 
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for  A® 

P  =J^{S,P„P^^  +  T-~S  for  drainage  ^^5<ol 

Pg-Pi  =Jo"^(^S,P,,P^'j+T-~S  for  imbibition  ^^5>ol 
from  which  one  can  conclude: 

K  (5,^,^)  =  K,  (5);  Kl  (S,P;,P;)  =  (1  - S) 

and  for>^^' 
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5  5  fd-  \ 

8S  ot  yot  j 

and 

drHs]  s  (8-  \ 

nf  -  d  = - =^-^s,  +T  — s.  for  imbibition  — 5  >  0 

8S  dt  '  [8t 

As  it  is  clear  from  the  expressions  for  forces  F,'^/  their  time  and  space  averages  are  zero  in  the 
first  order  in  A'.  It  follows  then  that  in  order  {5,(x,t),p;(x,r),pf  (x,r)}have  zero  time  and 
space  averages  the  following  choice  must  be  made; 

k;(5,:^,:^)=o;  Kf(s,:^,:^)=o 
:p((s,:^,:^)=o;  :pf(5,:^,:^)=o 
j;'”*(5)  =  0;  Jf{s)  =  0 

The  following  set  of  linear  equations  describes  then  small  vibrations  of  pressures  and  saturations 
in  the  first  order  in  A' : 

For  >1’ 


61 
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The  same  way  as  for  A\  ls,(s,F„P^,x,t),p'(s,P„P^,x,t),p^(s,P„F^,x,t)j  should  not  possess 
slow  non-periodic  component.  To  satisfy  that,  one  requests: 
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yf*  (5)  =  for  Imbibition  >  0 j 

Equations  for  {52(5>^.^.x,t),p^S,:^,:^,xy),p,^(5,:^,:^,xy))are  fully  similar  to  the  case 
for  A'  and  look  as  follows: 
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Here 

{.:.)={  )-pj={ )- 

Summarizing,  the  {s(x,r),^(x,t),Pg  (x,t)}  that  describe  slow  non-periodic  transport  in  second 
order  in  vibration  amplitude  satisly  the  following  set  of  equations: 
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Here  =  vf ^  |  and  7-^ V  -A=  =  0  as  5,  being  periodic. 
O  \  »3  y  1“|3  \\’-SJ 


a  (  5  1  K. 


Ar'-v(K,5'v(5P;)) 


dS  S 


s  ds 


;/;'v. 


^  as  ^  i-s  J(i-5)  '  ^  as  i-s  J''  ^■''s 


sk/i-s) -  K^i-s) - 


^  «'  V  -K,(S)  .(5)-'v(5.  .F,-)]-f!^,(S) '[v(s..').l(5.1.«) 


//:'-v. 


-k.(,-5) ^  f.«,(,-s)-V((,-5).'f/).(,-5)-'|((,-5).i^) 


aK^i-s)- 


(i-sr[v((i-s).f,')+l((,-5).,r;‘)’ 


P,-f<=y‘*(s)  +  ^  +T~S  for  drainage  f-S<0 

^  '  oS  2  dt  \dt 


Pg-Pi  =J"”  (5')  + - = - —  +  T — S  for  imbibition  ( — S>0 

^  ^  dS  2  dt  [dt 
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+s-b-s^ 
-e-b-s^  -V-U; 


^•V  -  +  r-V  -  'u,  +  f  -  'u 

f-V  -  'Uj  +  £--i-(l-5'j-V-  'U;  +  5-V- 


g 


As  it  has  been  mentioned  earlier,  ^F,Y  requires  next  order  correction  for  {u^,U;,Ug}.  They  are 
not  appear  in  the  equations  for  slow  variables  in  the  second  order  in  A  and  therefore  are  not 
computed  here. 
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Appendix  II. 

Approximate  solution  by  method  of  harmonic  balance. 

Assume  now  that  is  a  harmonic  wave: 


uj  ru 


U/  i  =  j  *u,  l-cos(ft)-t  +  k-x)  +  j  ‘u,  [•sin(ry-/  +  k-x) 


“.J  I 


and  calculate  additional  pressure  and  permeability  that  this  simple  vibration  induces  in  porous 
material.  First,  a  response  to  this  driving  force  needs  to  be  calculated: 


"p\  ^•cos(ry-/  +  k-x)  +  j  [•sin(iy-/  +  k -x) 


Substituting  this  into  equations  for  {5,(x,t),p,'(x,/),/?f  (x,t)}one  finds  the  following  linear 
system  for  amplitudes 


Waves  and  structural  vibrations. 

Keeping  the  terms  with  the  highest  order  in  {ry,k}  one  gets 


A 

S'' -P, 

0 

1 

0 

0 

-Tp, 

-1 

0 

. 

0 

. 

0 

0 

0 

0 

....> 

1 

0 

0 

-1 

. 

0 

-djjdS 

T-o) 

-1 

0 

1 

0 

-Tco 

-djjdS 

0 

-1 

0 

1 

where 
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k  JdK,  _ « 


k-K,-k  [dS  S  J  S  KdS  Sj  J  k-K,-k 


.  K)  VS-  fdK,  K, 


k-K-k  [  dS  l-Sjl-S"  KdS  l-S 


+  ^  •— P„-  -VP  +-^ 


iO 


K,  k  Kg  k 


q-k-‘'u^+r-k-‘n,+s-b-S-k-‘'Ug 
qk"u^  +  rk'}i,+sbSk'’Ug 
r  •  k  •  “"u,  +  £■  •  6  •  (l  -  5")  •  k  •  ^u,  +  5  •  k  •  ‘'Ug 
/•k-"u,+£--6-(l-5)-k-'u, +5-k-*Ug 


k  K.  k 


•k-K,->+L-y9,  +  /?g-(5''-l)-;0g,  •k.K,-(X-^u 


k  Kg  k 


k  Kg  k 


/?,.k.Kg.X+  (l-4 


•k-Kg-X+  A.-y^g+A-  (l--^)'  -1  Mg/  •k-Kg-(X-*»* 


Eliminating  \^p[, ' p[,  /pf  }  from  the  equations  above  one  gets  the  system  for  only: 


a/  ,  g  P,_ 

dS  S  1-5 


—T  •  G)+ Di  —  D 


T-(o-D,  +  D 


dJ_  ,  g  P,_ 

dS  S  1-5 


^  f'  fs 

^  fS  f' 


and5^=(^5,)  +(^5.)  =(^F'-^F^)  +(^F^-^F') 


Shaking. 

Spatial  wave  dependencies  are  neglected  in  this  case.  Keeping  only  the  terms  with  the  highest 
order  in  o>  one  gets: 
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Ot  '  '  } 


p,-(o  - 


v-K,-  *U;+  A  +  y^^/-— -(5’'- 1) 


(  vs 

V  K,— 

s 


p.  o}  ■■ 


V-K,- 


P,  (o  - 


v-K,-  X+  Ps^Psrf  \(^-s)  -A  -(X-X) 

I  V  ^  V 


Eliminating  p[, ' p[, /pf  }  from  the  equations  above  one  gets  the  system  for  only: 


dS  ^ 


-0)-\  T  +  - 


v(k,(i-s)-'vs) 

^0  1 
v(K,S’’VS) 


0)-\  T  +  - 


v(k,(i-s)‘’vs) 

1 

v(K,S'Vs) 


_^+p>  +p> 

dS  " 


F' 

^F‘  F^ 


[^F'  F^f  +(^ F^  +"  F'f 


ds  '  y  X  v(k,(i-5)"'vs  v(k,s-'vs) 
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Correlations  of  induced  saturation  and  pressure. 


Rewriting  equations  for  ‘p[, )  in  matrix  notations, 

L-’-s, +j-p;  =f; 

G".s,+j.pf  =f; 

C-s,-I-p{+I-pf  =0 
for  waves  and 
u^-s, +i.p;  =if 
G^*.s,-I.pf=f; 

C-s.-I-pJ+I-pf  =0 


for  shaking. 
Here: 


with  the  obvious  property  •  J  =  I; 


N, 

D, 


^.T\c 

/)J  ’  -{-T-co 


T  o)  ^ 
-dJjdSj 


From  this,  first,  expressions  for  |s,,p(,pf  jean  be  found: 

s,  =  (c + j .  L"  -  j  •  G"  )■'  (f;  -  f; ) 

pf  •(c+j.L--j.G')-‘(f-f;))  , 

for  waves  and 
s,  =(c+u''+G’*)'' 

p(  =  f/"  -  u"  •  (c + If  +  G )■'  (f;''  +  ) 

pf  =  -f  +  G  •  (c + u'’ + G  )■'  [i;" + rj’ ) 
for  shaking. 


And  then  the  correlations  entering  equations  for  the  second  order  might  be  expressed: 


2  T 

5,  =  s, -s,  = 


det(c  +  J-L’'’-J-G’") 

\.p;  = '(f; -f;)  "(c+j  L'-j  G-)'' .j  (f; -L- (c+j  L- -J  G')"' (r; -f;)) 

=  '(f;-f;)/(c+j-L--jG')''j.(f,'-G'.(c+j.L'-j.G-)''.(f;-f;)) 

^  (f;  -  f; )  /  (c  +  J  •  L”’  -  J  •  G”'  )■*  •  (f;  -  L”'  •  (c  +  J  •  L"  -  J  •  G"  )■'  •  (f;  -  f;  ))  •  k 
S^-Wpf  =V5,  -pf  =  "'sj-J-pf  •k  = 

"(f;-f;)-"(c+j-L’-'-j-G’-)"‘-(f;'-G'‘'-(c+j-L’*'-j-G'^)'‘-(f;'-f;))-k 


for  waves  and 


2  T 

S,  =  Si-S,  = 


det(c  +  L'''+G"") 

^  =  "s,  •  pj  ="  (f;  +  f; )  •"  (c + + G^'’  )■’  •  (f/"  -  u"  •  (c + + G^"  )■'  (f;  + 
^;t^= "s,  -pf  ="  (f;  +f;)-"  (c+L^'’  +G^")''  -(f;"  -g^^  -(c+u"  +G^'')‘‘(f;  +f;)) 


for  shaking. 
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Appendix  III 
Dynamics  of  saturation 

Detailed  sequence  of  snapshots  demontsrating  the  change  of  watert  saturation  with  time  is  shown 
below. 


Gas  and  liquid  pressure. 


Comparison  between  gas  and  liquid  pressure  with  structural  vibrations  is  shown  in  figures 
below.  One  can  see  the  degree  of  in-phase  oscillations. 
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T«i»r«6<»»  Snifter  H  »»)  OttpiKunw):  QwiH 


8 


Solid  versus  liquid  displacement:  the  mode  collinear  with  the  excitation. 

Next  set  of  figures  illustrates  the  comparison  between  displacement  of  solid  and  liquid  phases 
y-direction,  along  the  plate  thickness. 


